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Abstract  We  develop  an  ELL  AM  (Eulerian-Lagrangian  localized  adjoint  method)  scheme  to  solve  two-dimensional 
advection-dispcrsion  equations  with  all  combinations  of  inflow  and  outflow  Dirichlet ,  Neumann,  and,  flux  boundary 
conditions.  The  ELLAM  formalism  provides  a  systematic  framework  for  implementation  of  general  boundary  condi¬ 
tions,  leading  to  mass-conservative  numerical  schemes.  The  computational  advantages  of  the  ELLAM  approximation 
have  been  demonstrated  for  a  number  of  one- dimensional  transport  systems;  practical  implementations  of  ELLAM 
schemes  in  multiple  spatial  dimensions  that  require  careful  algorithm,  development  arc  discussed,  in  detail,  in  this  pa¬ 
per.  Extensive  numerical  results  are  presented  to  compare  the  ELLAM  scheme,  with,  many  widely  used  numerical 
methods  and  to  demonstrate  the  strength,  of  the  ELLAM  scheme. 


1  Introduction 

Many  difficult  problems  arise  in  the  numerical  simulation  of  advection-dispersion  equations, 
which  describe  the  transport  of  solutes  in  groundwater  and  surface  water,  the  displace¬ 
ment  of  oil  by  fluid  injection  in  oil  recovery,  the  movement  of  aerosols  and  trace  gases  in 
the  atmosphere,  and  miscible  fluid  flow  processes  in  many  other  applications.  In  indus¬ 
trial  applications,  these  equations  are  commonly  discretized  via  finite  difference  methods 
(FDM)  or  finite  element  methods  (FEM)  in  large-scale  simulators.  Because  of  the  enormous 
size  of  many  field-scale  applications,  large  grid-spacings  must  be  used  in  the  simulations. 
When  physical  dispersion  dominates  the  transport  process,  these  methods  perform  fairly 
well.  However,  when  advection  dominates  the  transport  process,  these  methods  suffer  from 
serious  numerical  difficulties.  Centered  FDM  (in  space  or  time)  and  corresponding  FEM 
often  yield  numerical  solutions  with  excessive  oscillations.  The  classical  space-upwinded  (or 
backward-in-time)  schemes  can  greatly  suppress  the  oscillations,  but  they  tend  to  generate 
numerical  solutions  with  severe  damping  or  a  combination  of  both.  Recent  developments  in 
effectively  solving  advection-dispersion  equations  have  generally  been  along  one  of  the  two 
approaches:  Eulerian  or  characteristic  methods.  Eulerian  methods  use  a  fixed  spatial  grid 
such  as  optimal  test  function  methods  of  Christie  and  coworkers  (1976),  Barrett  and  Morton 
(1984),  Celia  et  al.  (1989),  and  Bank  et  al.  (1990).  These  methods  attempt  to  minimize 
the  error  in  approximating  spatial  derivatives  and  yield  an  upstream  bias  in  the  resulting 
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numerical  schemes.  Hence,  they  are  susceptible  to  time  truncation  errors  that  introduce  nu¬ 
merical  dispersion  and  the  restrictions  on  the  size  of  the  Courant  number,  and  they  tend  to 
be  ineffective  for  transient  advection-dominated  problems.  They  generally  require  small  time 
steps  for  reasons  of  accuracy,  because  the  time  truncation  error  depends  on  liigh-order  time 
derivatives  of  the  solutions  that  are  large  when  a  sharp  front  passes  by.  Other  Eulerian  meth¬ 
ods,  such  as  the  Petrov-Galerkin  FEM  of  Westerink  and  Shea  (1989),  Bouloutas  and  Celia 
(1991),  and  the  total  variation  diminishing  scheme  of  Cox  and  Nishikawa  (1991),  attempt  to 
reduce  the  overall  truncation  error  by  using  negative  temporal  numerical  dispersion  to  cancel 
positive  spatial  numerical  dispersion.  Therefore,  they  also  suffer  from  the  Courant  number 
restrictions.  Also,  included  in  the  class  of  Eulerian  methods  is  the  streamline  diffusion  finite 
element  method  (SDM)  of  Hughes,  Brooks,  and  Johnson,  et  al.  (Hughes  and  Brooks  1979, 
Brooks  and  Hughes  1982,  Hughes  and  Mallet  1986,  Hughes  1995,  Hansbo  and  Szepessy  1990, 
Hansbo  1992,  Johnson  and  Navert  1981,  Johnson  1987,  Johnson  et  al.  1990,  Johnson  and 
Szepessy  1995,  Zhou  92,  and  Zhou  95).  Via  a  framework  of  space-time  FEM,  the  SDM  uses 
piecewise  polynomial  trial/test  functions  over  a  partition  on  a  space-time  domain  (spatial 
domain  x  current  time  interval).  By  defining  the  test  functions  delicately,  the  SDM  adds 
a  numerical  dispersion  only  in  the  direction  of  characteristics  (streamline)  to  suppress  the 
oscillation  and  does  not  introduce  any  cross-wind  diffusion.  Therefore,  this  method  possesses 
many  physical  and  numerical  advantages  other  Eulerian  methods  do  not  have.  However,  this 
method  contains  a  undetermined  parameter  in  the  test  functions  that  needs  to  be  chosen 
very  carefully  to  obtain  accurate  numerical  results.  If  the  parameter  is  chosen  too  small,  the 
numerical  solutions  will  exhibit  oscillations.  But  if  it  is  too  large,  the  SDM  will  introduce 
excessive  numerical  dispersion  and  seriously  smear  the  numerical  solutions.  Unfortunately, 
an  optimal  choice  of  the  parameter  is  not  clear  and  is  heavily  problem-dependent.  Moreover, 
the  number  of  unknowns  are  doubled  compared  to  many  standard  Eulerian  or  characteristic 
methods. 

Because  of  the  hyperbolic  nature  of  advective  transport,  characteristic  analysis  is  nat¬ 
ural  to  aid  in  the  solution  of  advection-dispersion  equations  and  has  led  to  many  related 
approximation  techniques,  including  the  method  of  characteristics  of  Carder  et  al.  (1964), 
Pinder  and  Cooper  (1970),  Benque  and  Ronat  (1982),  and  Hervouet  (1986);  the  charac¬ 
teristic  Galerkin  method  of  Varoglu  and  Finn  (1980);  the  Eulerian-Lagrangian  method  of 
Neuman  (1981);  the  transport-diffusion  method  of  Pironneau  (1982);  the  modified  method  of 
characteristics  of  Douglas  and  Russell  (1982),  and  Ewing  et  al.  (1983);  the  operator-splitting 
method  of  Espedal  and  Ewing  (1987),  Wheeler  and  Dawson  (1988),  and  Dalile  et  al.  (1990); 
and  the  Lagrangian-Galerkin  method  of  Morton  and  coworkers  (1988).  Characteristic  meth¬ 
ods  effectively  solve  the  advective  component  by  a  characteristic  tracking  algorithm  and  treat 
the  diffusive  term  separately.  These  methods  have  significantly  reduced  the  time  truncation 
errors  in  the  Eulerian  methods,  have  generated  accurate  numerical  solutions  even  if  large 
time  steps  are  used,  and  have  eased  the  Courant  number  restrictions  of  Eulerian  methods. 
Problems  with  many  characteristic  methods  arise  in  the  areas  of  rigorously  treating  bound¬ 
ary  fluxes  when  characteristics  intersect  inflow  or  outflow  boundaries  and  of  maintaining 
mass  conservation. 

The  Eulerian-Lagrangian  localized  adjoint  method  (ELLAM)  was  first  introduced  by  Celia 
et  al.  (1990),  Russell  (1990),  and  Herrera  et  al.  (1993)  for  the  solution  of  one-dimensional 
(constant-coefficient)  advection-diffusion  equations.  The  ELLAM  formalism  provides  a  gen- 
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eral  characteristic  solution  procedure  for  advection-dominated  problems,  and  it  presents  a 
consistent  framework  for  treating  general  boundary  conditions  and  maintaining  mass  con¬ 
servation.  Subsequently,  Russell  and  Trujillo  (1990),  Wang  (1992),  and  Wang  et  al.  (1992) 
derived  different  ELLAM  schemes  for  one-dimensional  linear  variable-coefficient  advection- 
diffusion  equations  with  general  inflow  and  outflow  boundary  conditions,  based  on  different 
(forward  or  backward)  techniques  for  the  tracking  of  characteristics  of  the  velocity  field.  Celia 
and  Ferrand  (1993),  and  Healy  and  Russell  (1993)  extended  ELLAM  to  a  finite-volume 
setting  for  one-dimensional  advection-diffusion  equations.  Ewing  (1991)  and  Dahle  et  al. 
(1995)  addressed  the  ELLAM  techniques  for  one-dimensional  nonlinear  advection-diffusion 
equations.  Ewing  and  Wang  (1991,  1993a,  1993b)  also  developed  ELLAM  schemes  for  the 
solution  of  one-dimensional  advection-reaction  equations  with  an  initial  condition  and  in¬ 
flow  boundary  condition.  In  addition,  Celia  and  Zisman  (1990),  and  Ewing  and  Wang 
(1993c,  1994)  generalized  ELLAM  schemes  for  one-dimensional  advection-diffusion-reaction 
transport  equations.  Wang  et  al.  (1995a),  and  Vag  and  coworkers  applied  ELLAM  schemes 
to  solve  the  systems  of  one-dimensional  reactive  transport  problems  from  bioremediation 
and  other  applications  (Vag  et  al.). 

While  the  computational  advantages  of  ELLAM  approximations  have  been  demonstrated 
for  one-dimensional  advection-dominated  problems  by  the  extensive  research  mentioned 
above,  practical  implementation  of  ELLAM  schemes  in  multiple  spatial  dimensions  requires 
careful  algorithm  development,  in  which  some  research  has  been  carried  out  in  this  direction. 
Russell  and  Trujillo  (1990)  addressed  various  issues  in  multidimensional  ELLAM  schemes. 
Wang  (1992)  developed  an  ELLAM  simulator  to  solve  two-dimensional  linear  advection- 
dispersion  equations  with  general  inflow  and  outflow  boundary  conditions  by  combining 
forward  and  backward  tracking  algorithms.  Theoretically  optimal-order  error  estimates  for 
the  derived  scheme  were  also  proved,  and  various  numerical  experiments  were  performed. 
Some  of  these  results  were  reported  in  Ewing  and  Wang  (1993a,  1993b)  and  in  Wang  et  al. 
(1996).  By  using  an  explicit  mapping  of  the  finite  elements  at  the  current  time  level  to 
the  spatial  grids  at  the  previous  time,  Binning  and  Celia  (Binning  1994,  Binning  and  Celia 
1996)  reported  on  a  finite-volume  ELLAM  formulation  for  unsaturated  transport  in  two 
dimensions.  Relations  and  differences  between  the  two  approaches  are  discussed  in  some 
detail  in  Section  4  of  this  paper.  In  an  upcoming  paper  Healy  and  Russell  also  developed 
a  finite-volume  ELLAM  scheme  for  two-dimensional  linear  advection-dispersion  equations 
(Healy  and  Russell).  Celia  (1994)  also  explored  the  development  of  an  ELLAM  scheme  for 
three-dimensional  advection-dispersion  equations. 

A  different  but  related  method  is  the  ‘characteristic-mixed  finite-element’  method  of  Ar- 
bogast  et  al.  (1992),  Yang  (1992),  and  Arbogast  and  Wheeler  (1995),  which  uses  piecewise- 
constant  space-time  test  functions.  As  with  the  standard  mixed  method,  a  coupled  system 
results  for  both  the  concentration  and  the  diffusive  flux.  The  theoretically  proven  error 
estimate  is  (0(Ax)3/"2)  for  grid  size  Ax,  which  is  suboptimal  by  a  factor  0((Arr)1,/2).  For 
ELLAM  schemes  with  piecewise  linear  trial/test  functions  for  one-dimensional  advection- 
diffusion  equations,  advection-diffusion-reaction  equations,  and  first-order  advection-reaction 
equations,  optimal-order  error  estimates  of  0(( Ax)2)  have  been  proven  by  Ewing  and  Wang 
(1991, 1993a, 1996),  Wang  et  al.  (1995a),  and  Wang  and  Ewing  (1995). 

Based  on  the  approach  presented  in  Wang  (1992),  an  ELLAM  scheme  is  developed  in  this 
paper  for  the  numerical  solution  of  two-dimensional  linear  advection-dispersion  equations 
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with  general  inflow  and  outflow  boundary  conditions.  We  have  organized  this  paper  as  fol¬ 
lows:  we  begin  in  Section  2  by  presenting  a  space-time  variational  formulation  of  the  model 
equations.  In  Section  3  we  derive  a  corresponding  ELLAM  scheme  for  this  formulation  with 
implementational  issues  discussed  in  Section  4.  Section  5  provides  a  brief  description  of 
some  well  studied  and  widely  used  methods,  including  the  Galerkin  finite  element  method, 
the  quadratic  Petrov-Galerkin  method,  the  cubic  Petrov-Galerkin  method,  and  the  stream¬ 
line  diffusion  finite  element  method.  In  Section  6  we  carry  out  numerical  experiments  and 
compare  the  performance  of  the  ELLAM  scheme  with  the  numerical  methods  described  in 
Section  5.  In  Section  7  we  summarize  our  observations  and  results. 

2  Variational  Formulation 

A  general  linear,  variable-coefficient  advection-dispersion  partial  differential  equation  in  two 
dimensions  can  be  written  as  follows: 

£tt(x,  t)  =  (i?(x,  t)u)t  +  V  •  ( V (x,  t)u  -  D(x,  t) Vm)  =  /(x,  t), 

(x,  t)  =  (x,  y,  t)  E  S  =  Q  x  J, 

where  ut  =  Vu  =  (|^,  j^)T ,  Q  is  a  spatial  domain,  J  =  [0,T]  is  a  time  interval.  The 
nomenclature  is  such  that  R(x.,t)  is  a  retardation  coefficient,  V(x,  t)  =  (Vi(x,  t),  V^x,  t)) 
is  a  fluid  velocity  field,  D(x,  t)  =  (D,j(x,i))T=1  is  a  diffusion-dispersion  tensor,  /(x,  t)  is 
a  given  forcing  function,  and  u(x,  t)  is  the  solute  concentration  of  a  dissolved  substance. 
Mathematically,  R  has  positive  lower  and  upper  bounds,  D(x,i)  is  a  symmetric  and  positive 
definite  matrix  with  uniform  lower  and  upper  bounds  that  are  independent  of  (x,  t). 

Let  the  space-time  boundary  T  =  dQ,  x  J  be  decomposed  as  the  union  of  an  inflow 
boundary  T*7*,  an  outflow  boundary  r(0),  and  a  noflow  boundary  T(  A)  (i.e. ,  T  =  P( ur(0)  U 
r(A)).  In  general,  an  inflow  boundary  during  one  time  period  might  become  an  outflow  or 
a  noflow  boundary  in  the  next  time  period  or  verse  versa.  At  P( or  r(0),  one  of  Dirichlet, 
Neumann,  or  Robin  (flux)  boundary  conditions  may  be  imposed  by  setting,  respectively, 

u(*,t)  =t/i')(x,t),  (x,<)er(i), 

— DVtt(x, t)  ■  n  =g^(x,t%  (x,i)er(,), 

(Vu  -  D Vii)(x,f)  •  n  =  /yV'lx./),  (x,i)  E  r(,), 

where  n  =  n(x)  is  the  outward  unit  normal,  i  =  /  or  O  represents  the  inflow  or  outflow 
boundary,  respectively.  A  noflow  boundary  condition  is  specified  at  T(A)  by 

(Vu  -  D Vii)(x,f)  •  n  =  0,  (x,i)  E  r(  V).  (2.3) 

In  addition  to  the  boundary  conditions,  an  initial  condition  u(x,  0)  =  u o(x)  is  needed  to 
close  Equation  (2.1). 

The  ELLAM  formalism  uses  a  time-marcliing  algorithm.  Let  Nt  be  a  positive  integer. 
We  define  a  partition  of  time  interval  J  =  [0,T]  by 

0  =  tg  <  t\  <  ^2  <  •  •  •  <  tn  <  .  .  .  <  1  <  ijVj  =  T. 


(2.2) 
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With  space-time  test  functions  w  that  vanish  outside  S„  =  Q  x  J„  with  J„  =  (f„_i,f„]  and 
are  discontinuous  in  time  at  time  | .  one  can  write  a  space-time  variational  formulation 
for  Equation  (2.1)  as  follows: 


u  w)(x,tn)  dx  +  Vw  ■  (DVti)  dxdt 
+  f  (V u  —  DYn)  -n  w  dS  —  (  u  (Rwt  +  V  ■  Vw)  dx.dt 
=  fa  ( R  u  tc)(x,(+_,)  <ix  +  /E  /  w  dxdt, 


(2.4) 


where  T„  =  dQ,  x  J„  and  r«(x,  =  lim  u;(x,  t). 


In  the  ELLAM  framework,  one  should  define  the  test  functions  w  to  satisfy  the  equation 
Rwt  +  V  •  Vw  =  0,  so  the  last  term  on  the  left-hand  side  of  Equation  (2.4)  vanishes.  However, 
in  general  one  cannot  track  characteristics  exactly  for  a  variable-velocity  field.  Nevertheless, 
this  adjoint  term  should  be  small  if  one  can  track  the  characteristics  reasonably  well,  and  the 
test  functions  are  constant  along  the  approximate  characteristics.  In  fact,  we  have  proved 
an  optimal-order  convergence  rate  for  the  derived  ELLAM  scheme  even  if  a  one-step  Euler 
algorithm  is  used  in  the  characteristic  tracking  and  this  adjoint  term  is  dropped  (Wang  1992, 
Ewing  and  Wang  1996).  To  conserve  mass,  the  test  functions  should  sum  to  one  (Celia  et  al. 
1990).  The  scheme  developed  in  this  paper  satisfies  this  condition.  In  this  case  dropping  the 
last  term  on  the  left-hand  side  does  not  affect  mass  conservation  since  that  term  vanishes  if 
w  =  1  (Russell  and  Trujillo  1990,  Wang  1992). 


We  are  now  in  a  position  to  rewrite  Equation  (2.4).  Given  a  point  (x,  1),  with  t  E  [£„_i,  tn\, 
we  consider  the  initial-value  problem  for  the  ordinary  differential  equation 


dx 

dt 

x(f) 


=  V*(x,f) 

=  x, 


V(x,t) 

R(x,t)  ’ 


(2.5) 


which  tracks  the  characteristics  from  (x,t).  We  denote  the  solution  of  this  equation  at  time 
6  E  Jn  by  X(0;x,F)  (Healy  and  Russell  1993).  This  notation  can  refer  to  tracking  either 
forward  or  backward  in  time.  In  particular,  we  define 


x*  =  X(t„_1;  x,  tn), 

V  ’  ’  h  (2.6) 

x  -  X(/,.:x.  C  i). 

Thus,  (x,i„)  backtracks  to  (x' .  t„  \ )  and  (x,  G-i)  tracks  forward  to  (x,t„)-  In  the  numerical 
scheme,  an  exact  tracking  is  preferred  whenever  possible.  However,  it  is  impractical  in  most 
applications.  In  practice,  one  can  use  either  a  one-point  Euler  quadrature,  a  multiple  micro- 
time-step  tracking  within  a  global  time  step,  or  a  Runge-Kutta  quadrature  in  the  tracking 
of  characteristics.  Note  that  in  many  applications,  Equation  (2.1)  is  usually  coupled  with  an 
associated  potential  or  pressure  equation,  whose  solution  is  often  obtained  via  the  mixed  finite 
element  method.  In  this  case,  a  Raviart-Thomas  space  is  often  used  for  the  velocity  field, 
which  is  calculated  at  each  cell  interface.  Within  each  cell,  Vj(x,£)  (or  V<i (x,  £))  is  piecewise 
linear  (or  constant)  in  the  x  direction  and  piecewise  constant  (or  linear)  in  the  y  direction. 
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Under  the  assumption  that  the  velocity  field  is  steady,  a  semi-analytical  technique  has  been 
developed  by  Pollock  (1988),  Schafer-Perini  and  Wilson  (1991),  and  Heally  and  Russell 
(1993),  to  track  the  characteristics  on  a  cell-by-cell  basis.  Recently,  Lu  (1994)  extended 
this  semi-analytical  approach  to  non-steady  velocity  fields  where  velocity  is  assumed  to  vary 
linearly  in  time  within  each  time  interval. 

To  accurately  measure  the  time  period  taken  for  a  particle  to  move  along  a  characteristic 
from  the  previous  time  level  (or  from  the  inflow  boundary)  to  the  current  time  level  (or 
the  outflow  boundary),  we  introduce  space-time  location-dependent  time  steps.  We  use 
Figure  1(a)  to  illustrate  how  these  are  defined.  In  this  figure,  we  use  letters  A-D  to  denote 
points  (i)  at  the  future  time  step  tn  (points  B  and  C)  or  (ii)  on  the  outflow  boundary 
(points  A  and  D),  while  A*-D*  will  denote  the  corresponding  feet  of  their  characteristics. 
In  our  example,  we  have  set  the  rear  planes  x  =  a  and  y  =  c  as  the  inflow  boundaries,  with 
the  frontmost  planes  x  =  b  and  y  =  d  as  the  outflow  boundaries.  Time  is  represented  in 
the  vertical  direction.  For  any  x  £  ft  at  time  t ,,  ,  we  define  a  time  step  Ad (x)  =  At  = 
tn  —  i  if  the  characteristic  X(0;x,  tn)  does  not  backtrack  to  the  space-time  boundary  F„ 
during  the  time  period  Jn  (this  case  is  illustrated  by  point  B  at  time  tn  in  Figure  1(a))  and 
A fd)(x)  =  tn  —  f*(x)  otherwise  (  see,  for  example,  point  C).  In  the  latter  case  where  the 
foot  of  the  characteristic  (point  C *)  lies  on  an  inflow  boundary,  f*(x)  G  Jn  is  the  time  when 
X(0;x,  £„)  intersects  the  boundary  Yn  (i.e,  X(£*(x);  x,£„)  G  T„).  Similarly,  for  any  point  on 
the  outflow  boundary  (x,  t)  E  F 1  (e.g.,  points  A  or  D),  we  define  A (x,  t)  =  t—tn_ i  if  the 
characteristic  X(0;  x,  t)  does  not  intersect  F„  during  the  time  period  t\,  otherwise  we  set 
A (x,  t)  =  t  —  f*(x,  t).  The  first  case  is  illustrated  by  point  A  on  the  space-time  boundary 
X  =  b.  while  the  second  case  is  demonstrated  by  point  D  on  the  space-time  boundary  y  =  d. 
Here  we  denote  by  t*(x,t)  G  the  time  when  X(0;x,  i)  intersects  T„. 

By  enforcing  the  backward  Euler  quadrature  at  the  current  time  tn  and  at  the  outflow 
space-time  boundary  F;)^ ,  we  approximate  the  space-time  volume  integral  of  the  source  term 
(the  second  term  on  the  right-hand  side)  in  Equation  (2.4)  by  an  integral  at  time  tn  and 
one  at  rj,0*  by  following  the  characteristics.  Here  F';.'1  =  fl  Jn  (■ i  =  I,0,N)  represents 
the  space-time  inflow,  outflow,  and  noflow  boundaries  during  the  time  interval  ./„ .  To  avoid 
confusion  in  the  following  derivation,  we  replace  the  dummy  variables  x  G  O  and  t  E  Jn 
in  this  term  by  y  £  ft  and  6  E  Jn.  Thus,  J'yj  f  w(x,  t)  dxdt  =  f  w(y,6)  dydO.  Let 

E 1  C  E„  be  the  set  of  points  in  the  space-time  strip  E„  that  will  flow  out  of  E„  during 
the  time  interval  Jn.  We  decompose  E„  to  be  the  union  of  E(°:i  and  E„  —  .  For  any 

(y,  6)  E  E,j  —  there  exists  a  point  x  £  ft  such  that  x  =  X ( t,„ :  y,  6).  Thus,  we  can  invert 
this  relation  to  obtain  y  =  X(0;x,  tn).  Similarly,  for  any  (y,0)  E  E 1 ,  there  exists  a  pair 
(x,t)  E  rf>  such  that  y  =  X(0;x,  t).  By  splitting  the  space-time  volume  integral  on  En  as 
one  on  E„  —  E),0*  and  one  on  E;,0)  and  applying  the  backward  Euler  quadrature  at  time  tn 


6 


for  the  first  and  at  boundary  V\(i  >]  for  the  second,  we  obtain  the  following  equation 


/( y,0)w(y,6)  dydd 

n 

=  J  (0)  f(y,d)w(y,d)  dyd6  +  j  (Q)  f(y,d)w(y,d)  dydd 
=  /s  s(OJ  f(X(d-x,tn),d)  w(X(8;x,tn),8)  dXdd  (2.7) 

+  J^0)f(X(d;x,t),d)  w(X(8]  x,  t),  d)  dXdd 
=  I^At{n(x.)  fn  wn(x)  dx  +  j^(Q)  A t{0)(x,t)  f  w  V  •  n  dS  +  Ef(w). 


Here  /„(x)  =  f(x,tn),  Ef  is  the  truncation  error  from  the  application  of  the  backward  Euler 
quadrature.  In  the  derivation  of  Equation  (2.7),  we  have  used  the  fact  that  the  test  function 
w  is  constant  along  the  characteristics. 


Similarly,  we  can  rewrite  the  diffusion-dispersion  term  as 


Yw  ■  (DYu)(y,d)  dydd 

n 

=  j  At(7)(x)  Ywn  ■  (D „Vu„)(x)  dx  (2.8) 

+  J  (0)  At{0)  (x,  t)  Yw  ■  (DVtt)  V  •  n  dS  4-  Ejy(u,  w ), 


where  Ejy(u,w)  is  the  truncation  error  term. 

Substituting  Equations  (2.7)  and  (2.8)  for  the  second  terms  on  both  the  left-  and  right- 
hand  sides  of  Equation  (2.4),  we  obtain  the  following  variational  formulation 


Wn  dx+  I^At(n{x)  Ywn  ■  (DnYun)  dx 
4-  f  {a)  A  t<°>  (x,  t)  Yw  ■  (DV«)  V  •  n  dS  +  f  (Yu-  DVti)  •  n 

~ku(R  wt  +  V  •  Yw)  dxdt 
=  j  Rn_iUn_i  w+_i  dx  +  j  At{E  (x)  /„  wn  dx 
-1-  jy(0)  At{0)(x,  t )  /  w  V  •  n  dS  4  E(u ,  w), 
where  E(u,w)  =  —Ej}(u,w)  4-  Ef(w). 


w  dS 


(2.9) 


3  An  ELLAM  Scheme 

While  the  numerical  scheme  can  be  derived  for  a  general  domain  Q  with  a  quasi-uniform 
triangular  or  rectangular  partition,  we  assume  the  domain  Q  =  ( a,b )  x  (c,d)  for  simplicity 
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to  be  a  rectangular  domain  with  a  uniform  rectangular  partition: 


Xj  =  a  +  iAx, 

Vj  =  c  +  jAy, 

where  Nx  and  Ny  are  two  positive  integers.  We  dehne  the  test  functions  wiy  to  be  piecewise- 
linear  “hat”  functions  at  time  t„  ( )  =  S^Sji  where  x/,/  =  (.r;..  //.•).  8ik  =  1  if  i  =  k 
and  0  otherwise)  and  to  be  constant  along  the  characteristics.  At  time  tn,  we  also  use 
piecewise-linear  trial  functions  U(x,tn). 


b  —  a 

i  =  0,1,.. 

. ,  Nx, 

Ax  = - 

Nx 

d  —  c 

J  =  0, 1, . . 

N 

> 

II 

F 

3.1  Interior  Nodes  and  Noflow  Boundary 

In  this  subsection  we  develop  the  scheme  at  the  nodes  inside  Q,  or  on  the  noflow  boundary 
that  are  neither  related  to  the  inflow  boundary  F';/ 1  nor  the  outflow  boundary  F 1 .  It 
is  assumed  that  the  type  of  boundary  (inflow,  outflow,  or  noflow)  will  be  kept  unchanged 
during  the  time  interval  ./„ .  Let 


r  n(q)  =< 

[(x,t)  G  r„  X  =  q ^ 

\  F< 

x,<)  x  =  q,y 

r n(q)  =  {(x,<)  ETny  =  q] 

define  the  C  our  ant  numbers 

•  F< 

x,i)  x  G  [a,  t 

Cu = 

max 

|Vi(x,<)|A< 

X 

(X 

•t)eF „  (q) 

Ax 

Cu\q]  = 

max 

\V2(x,t)\At 

A 

for  q  =  a,  b, 
for  q  =  c,  d. 


q  =  a,  6, 
q  =  c,  d. 


(3.1) 


(3.2) 


If  Tn(q)  (q  =  a,b )  is  an  inflow  or  outflow  boundary  (which  implies  that  Cv\q]  >  0),  we 
define  1C ^  to  be  IC\q'1  =  [Cu^]  +  1,  where  [Cu[^\  is  the  integer  part  of  Cu[qK  If  r„(g)  is 
a  noflow  boundary  (which  indicates  that  Cu\q *  =  0),  we  define  1 =  0.  1C ^  is  defined 
similarly.  In  addition,  we  define 


n  =  [a  +  IC^Ax ,  b  -  IC^Ax]  x  [c  +  IC^Ay,  d  -  IC^Ay], 

^ ij  i  ?  x  \yj— 1?  %+l]  ‘ 

Furthermore,  let  EF-  be  the  prism  obtained  by  backtracking  0,y  along  the  characteristics 
from  tn  to  F-i,  and  E „ rj  to  be  the  prism  obtained  by  tracing  Q (/  forward  along  the  charac¬ 
teristics  from  tn_ i  to  tn. 

When  f Ijj  C  D,  EF-  or  E„,;j  does  not  intersect  Tjd  or  rj,0*  during  the  time  period  Jn. 
The  third  and  fourth  terms  on  the  left-hand  side  and  the  third  term  on  the  right-hand  side 
of  Equation  (2.9)  vanish.  Dropping  the  last  terms  on  both  sides  of  Equation  (2.9),  replacing 


(3.3) 


the  exact  solution  u  and  the  general  test  function  w  by  the  piecewise-linear  trial  function  U 
and  test  function  we  obtain  the  following  equation 


jn  R„U„  Wijn (x)  dx+  J  At  Vwijn  •  (D „VU„)(x)  dx 
=  tr'i  h  ,(x)  dx+  At  /„  my„(x)  dx, 


with  iC;  :„(x)  =  U'ij(x,  t„).  Note  that  in  (3.4),  the  integrals  at  time  are  actually  defined 
on  Q (/  (with  the  obvious  modification  near  the  boundary  dQ)  since  Q rj  is  the  support  of 
m,?„.  The  first  term  on  the  right-hand  side  is  actually  defined  on  the  backtracked  image 
(at  time  t„_i)  of  Q (/  at  time  ,  which  can  be  of  a  very  complicated  shape  and  not  aligned 
with  any  elements  in  fi  at  time  tn_1  due  to  the  effect  of  the  velocity  field  even  though  fi , (  is 
rectangular.  Consequently,  the  evaluation  of  this  term  is  tricky  and,  in  fact,  crucial  to  the 
accuracy  and  mass  conservation  property  of  the  scheme.  This  will  be  discussed  in  detail  in 
Section  4.  At  this  point,  one  can  easily  see  that  the  scheme  has  a  9-banded,  symmetric  and 
positive-definite  coefficient  matrix. 

3.2  Inflow  Boundary  Conditions 

In  contrast  to  many  characteristic  methods  that  treat  boundary  conditions  in  an  ad  hoc 
manner,  the  ELLAM  scheme  naturally  incorporates  boundary  conditions  into  its  formulation. 
Thus,  one  can  approximate  boundary  conditions  accurately.  In  fact,  if  •  intersects  the 
inflow  boundary  T;/1 ,  the  test  function  w,j  assumes  nonzero  values  on  portions  of  T„ .  Thus, 
the  fourth  term  on  the  left-hand  side  of  Equation  (2.9)  does  not  vanish.  For  an  inflow  flux 
boundary  condition,  the  scheme  becomes 


l^RnU„  wijn (x)  dx  +  I^At(T)(x)  V W,:jn  ■  (D„VC/„)(x)  dx 

=  ,  (x )  dx+  I^At{T](x)  fn  w„(x)  dx  (3.5) 

-  (I)  g(3I]  w,u(x,t)  dS. 

Keep  in  mind  that  the  first  term  on  the  right-hand  side  of  Equation  (3.5)  is  now  defined 
on  the  image  (at  time  tn_1)  of  the  portion  of  Q,y  that  is  not  taken  to  the  boundary  T„. 
The  part  of  the  integral  that  is  missing  from  this  term  is  picked  up  by  the  last  term  on  the 
right-hand  side  of  Equation  (3.5),  which  is  defined  on  the  image  of  the  portion  of  Q rj  which 
is  taken  to  the  boundary  T„.  Notice  that  the  factor  A t^(x)  at  time  tn  now  depends  on  x, 
since  X(0;  x,  /„)  can  encounter  the  boundary  T„.  At1'1'1  (x)  reflects  the  time  period  over  which 
the  diffusion-dispersion  and  source  act.  For  an  inflow  flux  boundary  condition  the  derived 
scheme  still  has  a  9-banded,  symmetric  and  positive-definite  coefficient  matrix. 

Repeating  the  above  derivation  for  an  inflow  Diriclrlet  boundary  condition  yields  the 
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following  equation 


jQ  R  JJn  wijn(x)  dx  +  J  At(/r)(x)  Vwijn  ■  (D„W„)(x)  dx 

-  /(7)  (DVC/)  •  n  Wij(x,t)  dS 

=  ^  Rn-\Un-i  ^„_!(x)  dx  +  J^At{I)(x)  fn  Wjjn{x)  dx 

-  V  •  n  <yi/'  Wi-ix.l.)  dS. 

72 


(3.6) 


While  all  other  terms  in  Equation  (3.6)  are  similar  to  those  in  Equation  (3.5),  the  third 
term  on  the  left-hand  side  couples  the  unknown  boundary  diffusive  flux  with  unknown  in¬ 
terior  function  values.  If  one  simply  represents  VU  as  a  discrete  gradient  dependent  on 
imposed  boundary  values  of  U,  one  might  introduce  strong  temporal  truncation  errors.  To 
overcome  this  difficulty,  we  approximate  VC/(x,  t)  at  the  inflow  boundary  E;/1  implicitly  by 
V£/(X(f„ ;  x,  I: ) ,  t„ )  at  time  tn .  This  removes  the  difficulty  of  evaluating  an  unknown  diffusive 
boundary  flux.  The  error  introduced  is  small  since  it  is  along  the  characteristics  and,  in  fact, 
does  not  affect  the  convergence  rate  of  the  scheme  (Wang  et  al.  1995).  Note  that  this  term 
introduces  non-symmetry  to  the  coefficient  matrix  near  the  inflow  boundary. 

As  with  the  standard  finite  element  methods,  the  Diriclilet  boundary  condition  is  essential 
and  is  imposed  directly  on  the  solution  u  with  no  degrees  of  freedom  on  the  inflow  boundary 
rj/h  However,  the  test  functions  should  sum  to  one  to  conserve  mass  (Celia  et  al.  1990). 
Thus,  on  each  element  Q (/  =  rc,;+i]  x  [yj-i,  Vj+i]  that  has  at  least  one  vertex  on  the  inflow 

boundary  T;/1  the  test  functions  are  chosen  such  that  they  sum  to  one  on  this  element.  For 
example,  suppose  x  =  xq  =  a  is  an  inflow  boundary.  Then  the  interior  nodes  xij  =  (aq  ,yj) 
with  X\  =  a  +  Ax  and  1  <  j  <  Ny  —  1  are  adjacent  to  the  inflow  boundary  x  =  a.  In 
this  case,  the  corresponding  test  functions  must  satisfy  w1j(x,y)  =  W\ j(x\..  //).  i.e.  they  are 
constant  in  x  direction  over  the  interval  [a,aq]. 

A  similar  derivation  to  that  of  Equation  (3.5)  yields  a  scheme  for  Equation  (2.1)  with  an 
inflow  Neumann  boundary  condition.  This  differs  from  Equation  (3.5)  in  that  g':[]  is  replaced 
by  g.2^  and  an  extra  term  /yqq  U  wt  j  V  •  n  dS  appears  on  the  left-hand  side  of  the  equations. 

Because  V  •  n  Ip (7)  <  0,  this  term  has  a  different  sign  from  the  first  term  on  the  left-hand 

n 

side  of  Equation  (3.5). 

If  Tlr)  can  be  decomposed  as  Tj/*  =  T^  U  U  T^  where  inflow  Dirichlet,  Neumann, 

and  Robin  boundary  conditions  are  imposed  on  T*/1, ,  T*/),  and  T*) j,  respectively,  one  can 
write  out  the  scheme  accordingly. 

3.3  Outflow  Boundary  Conditions 

The  situation  at  the  outflow  boundary  V\(i  >]  is  different  from  that  at  an  inflow  boundary  Tj/*. 
The  number  of  spatial  degrees  of  freedom  crossing  the  outflow  boundary  T;,01  is  essentially  the 
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Courant  number  in  the  normal  direction.  To  preserve  the  information,  one  should  discretize 
in  time  at  the  outflow  boundary  T;,01  with  about  the  same  number  of  degrees  of  freedom. 
More  precisely,  we  define 


Cu (0)  =  max 

(X, fieri0) 


|Vi(x,<)|A<  |l^(x,  t)\At  1 


Ax 


Ay 


) 


(3.7) 


and  7(7*°)  =  [C'utUj]  +  1.  Then  we  define  a  uniform  local  refinement  in  time  at  the  outflow 
boundary  Tj0) 

iAt 

tn  i  tn  ~]~(j  for  %  0?  1?  •  •  •  ?  7(7. 


Of  course,  if  one  is  not  interested  in  accurate  simulation  near  the  outflow  boundary  Ti0), 
one  need  not  refine  in  time  at  T 1 .  This  corresponds  to  the  choice  of  7(7  =  1.  In  any  case, 
we  define  the  test  functions  to  be  the  piecewise-linear  “hat”  functions  at  the  nodes  at 
the  outflow  boundary  T 1  and  to  be  constant  along  the  characteristics.  We  define  the  trial 
functions  T/(x,  t)  for  (x,/)  E  V\(I>]  to  be  the  piecewise-linear  functions  at  V 1 .  Incorporating 
the  outflow  Neumann  boundary  condition  into  Equation  (2.9)  yields  a  scheme  for  Equation 
(2.1)  with  the  stated  boundary  condition  as  follows: 

Jq  RnUn  wijn (x)  dx  +  I^At{n(x)  Vwijn  ■  (D„V7/„)(x)  dx 

+  f  A/'"1  (x.  /)  V/r.,:  •  (DVT/)  V  •  n(x./)  dS  +  f  U  wt]  V  •  n(x./)  dS 

=  I  7?„_iT/„_i  tn'-,.  ,(x)  dx+  I  At{T](x)  fn  wijn(x)  dx 

+  JT(0)  A/'(,)(x.  /)  /  Wij  V  •  n(x, /)  dS  -  J^(Q]  gf]  dS. 


(3.8) 


Because  T/,  not  VC7,  is  defined  as  unknowns  at  the  outflow  boundary  T 1 ,  it  is  difficult  to 
approximate  VC7  •  nLoj  numerically.  To  circumvent  this  difficulty,  we  utilize  the  boundary 

n 

condition  Equation  (2.2)  to  express  VC7  •  n|^(o)  in  terms  of  T/L(oj  and  the  tangential  com- 

ponent  of  VT/|^(0),  which  can  be  computed  by  differentiating  U  on  T;,01 .  To  demonstrate 

the  ideas,  we  assume  that  T„(6)  (i.e,  the  “eastern”  boundary  x  =  b)  is  an  outflow  boundary. 
The  outflow  Neumann  boundary  condition  in  Equation  (2.2)  now  reads 


-DVm  •  n  =  7)  1 1 


D12uy  =  g\ 


( o ) 


1 


from  which  one  can  express  V«  •  n 
in  this  case)  and  u  as  follows 


■u 


X 


u:r  in  terms  of  the  tangential  component  of  Vtt  (uy 

P\2  uy  +  gj  1 


This  yields 


BVu  =  - 


( o ) 
92  , 


ID 


uy  -  D2ig2 


(o)  ’ 


D 


n 
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where  |D|  =  D11D22  —  D12D21  is  the  determinant  of  D. 

Using  the  facts  that  the  test  functions  Wjj  satisfy  the  equation  Rwt  +  V  •  Vw  =  0  and  that 
V’|  |-p(/,j  >  0,  since  h;)’1  is  an  outflow  boundary,  we  can  denote  Vuj;j  •  n|r(/,)  in  the  third  term 

on  the  left-hand  side  of  Equation  (3.8)  by  wx  =  —(Riut  4-  VoWy) /V\ .  Then  we  can  rewrite 
the  third  term  on  the  left-hand  side  of  Equation  (3.8)  as 


I  A *<°>(x,<)  VwtJ  •  (DW)  V  •  n(x,i)  dS 


=  JT(b]  Af(0)  (x,  t)  WijyUyi^t)  dS 


(x,  t)dS. 


(3.9) 


Substituting  this  equation  for  the  third  term  on  the  left-hand  side  of  Equation  (3.8),  we 
obtain  a  numerical  scheme  for  Equal  ion  (2.1)  with  an  outflow  Neumann  boundary  condition. 
The  derived  scheme  has  a  symmetric  and  positive-definite  coefficient  matrix. 


Since  the  numerical  solution  U  is  known  at  time  tn_1  from  the  computation  at  time  tn_1, 
there  are  no  degrees  of  freedom  on  the  boundary  T\(I  >]  at  time  To  conserve  mass,  the 

test  functions  on  T ^  that  intersect  Q  at  time  1  are  chosen  such  that  they  sum  to  one; 
this  was  discussed  following  Equation  (3.6). 


Incorporating  the  flux  boundary  condition  into  Equation  (2.9),  one  can  derive  a  scheme 
similar  to  Equation  (3.8)  except  that  the  last  term  on  its  left-hand  side  disappears  and  is 
replaced  by  .  Again,  we  need  to  express  Vu  ■  nLoj  and  Vm,-7  •  n|-p(0)  by  their  tangential 


derivatives  and  functional  values.  If  we  still  assume  that  T;)'1  is  an  outflow  boundary,  the 
outflow  flux  boundary  condition  in  Equation  (2.2)  now  becomes 


Dnux  +  D12uy  =  hi  u  -  gf\ 


which  yields 


Ur  = 


D 12  Vl  u  -  g3 

Uy  + 


(O) 


D 


11 


Combining  these  two  equations  gives 


]DI 

Du 


D 11 

D 


21 


D  V*  =  (  Vm  -  gf] ,  L±uy  +  -2-(VlU  -  gf]) 

V 11  3 


Similar  derivation  to  Equation  (3.9)  results  in  the  following  equation 
JT{b)  At(°>(x,t)  V uijj  ■  (D VU)  V  •  n(x,  t)  dS 


Lb)  At{0)(x,t) 

n 

+ h 


hq  D 


wijyUy  ~b 


V?D- 


21 


Du 

(RW,jjt  +  VgWjjy)  ~ 


-w 


Du 

V1D21W 


vi  y 


D 


11 


U  —  hi  (Rwjjt  +  hyu.)  ijy)U 
(x,  t)dS. 


(x,  t)dS 


(3.10) 
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Substituting  Equation  (3.10)  for  the  third  term  on  the  left-hand  side  of  Equation  (3.8), 
dropping  the  last  term  on  its  left-hand  side,  and  replacing  by  we  obtain  a  numerical 
scheme  for  Equation  (2.1)  with  an  outflow  flux  boundary  condition. 

For  Equation  (2.1)  with  an  outflow  Dirichlet  boundary  condition,  the  equations  at  the 
outflow  boundary  define  the  unknowns  to  be  the  normal  derivatives  of  the  solutions  and 
are  decoupled  from  the  equations  at  the  interior  domain  given  by  Equation  (3.4).  They  are 
omitted  here  since  they  are  only  needed  for  mass  conservation. 

If  E 1  can  be  decomposed  as  Rj/*  =  V U  T U  where  outflow  Dirichlet,  Neumann, 

and  Robin  boundary  conditions  are  imposed  on  T^,  ru)! .  and  T ,  respectively,  one  can 
write  out  the  scheme  accordingly. 

4  Implementation 

4.1  Evaluation  Of  Integrals  and  Tracking  Algorithms 

Some  integrals  in  the  numerical  scheme  derived  in  Section  3  are  standard  in  FEM  and  can 
be  evaluated  fairly  easily,  while  others  can  be  difficult.  In  this  subsection,  we  discuss  the 
evaluation  of  the  integrals  in  Equation  (3.4)  and  will  discuss  the  treatment  of  boundary 
terms  in  Equations  (3.5)-(3.10)  in  the  next  subsection. 

Note  that  the  trial  function  C(x,  tn)  and  test  functions  m;j(x,  tn)  are  defined  as  standard 
tensor  products  of  piecewise-linear  functions  at  time  the  integrals  in  Equation  (3.4)  are 
standard  in  finite  element  methods  except  for  the  first  term  on  the  right-hand  side.  In  this 
term,  the  value  of  [/(x,i„_i)  is  known  from  the  solution  at  time  tn_  1.  However,  keep  in  mind 
that  the  test  functions  w^n_1  =  limi_j>i+  w =  /r.;(x.  /  ,>••).  where  x  =  X(R;  x,  R-i)  is 
the  point  at  the  head  corresponding  to  x  at  the  foot.  The  evaluation  of  this  term  becomes 
much  more  challenging  in  multiple  dimensions,  due  to  the  multi-dimensional  deformation  of 
each  finite  element  Q (/  on  which  the  test  functions  are  defined  as  the  geometry  is  backtracked 
from  time  tn  to  time  tn_  1. 

In  modified  method  of  characteristics  and  many  other  characteristic  schemes,  this  term 
has  traditionally  been  rewritten  as  an  integral  at  time  tn,  with  the  standard  value  of  •?/;,, (x,  tn) 
but  backtracking  to  evaluate  £/(x*,i„_i)  where  x*  =  X(£„_i;x, tn)  is  the  point  at  the  foot 
corresponding  to  x  at  the  head  (Dahle  et  al.  1990,  Douglas  and  Russell  1982,  and  Espedal 
and  Ewing  1987).  As  a  matter  of  fact,  it  has  been  shown  that  in  characteristic  methods  the 
backward  tracking  algorithm  is  critical  in  the  evaluation  of  this  term,  which  is  in  turn  crit¬ 
ical  to  the  accuracy  of  the  scheme  (Baptista  1987).  Because  of  this,  the  backward  tracking 
algorithm  has  been  used  in  many  ELLAM  works  (Binning  1994,  Celia  et  al.  1990,  Dahle  et 
al.  1995,  Ewing  and  Wang  1991,  1993a,  1993b,  1994,  Russell  1990,  Wang  et  al.  1992,  1995). 
However,  for  multidimensional  problems  the  evaluation  of  this  term  with  a  backtracking 
algorithm  requires  significant  effort,  due  to  the  need  to  define  the  geometry  at  time  i„-i, 
which  requires  mapping  of  points  along  the  boundary  of  the  element  and  subsequent  inter¬ 
polation  and  mapping  onto  the  fixed  spatial  grid  at  the  previous  time  level  tn_  1.  Binning 
and  Celia  (Binning  1994,  Binning  and  Celia  1996)  used  such  a  mapping  in  two  dimensions 
in  a  procedure  that  was  computationally  very  intensive,  especially  when  part  or  all  of  the 
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element  being  mapped  intersects  a  space-time  boundary  F„ .  This  approach  is  considered 
impractical  in  two  and  three  dimensions  (Binning  1994,  Celia  1994).  For  one-dimensional 
problems,  the  evaluation  of  this  term  is  relatively  simple  since  the  boundaries  of  the  spatial 
elements  are  points  rather  than  lines  or  surfaces.  In  this  case,  these  problems  were  overcome 
in  the  works  cited  above. 

The  most  practical  approach  for  evaluating  this  term  is  to  use  a  forward  tracking  algo¬ 
rithm,  which  was  proposed  by  Russell  and  Trujillo  (1990),  and  was  implemented  by  Heally 
and  Russell  for  a  one-dimensional  problem  (Healy  and  Russell  1993),  and  by  Ewing  and 
Wang  (1993a, 1993b),  and  Wang  (1992)  for  a  two-dimensional  problem.  This  would  enforce 
the  integration  quadrature  at  tn-\  with  respect  to  a  fixed  spatial  grid  on  which  Rn-\  and 
[/„_ i  are  defined,  the  difficult  evaluation  is  the  test  function  wf-  n_1.  Rather  than  backtrack¬ 
ing  the  geometry  and  estimating  the  test  functions  by  mapping  the  deformed  geometry  onto 
the  fixed  grid,  discrete  quadrature  points  chosen  on  the  fixed  grid  at  t„--\  in  a  regular  fashion 
(say,  standard  Gaussian  points)  can  be  forward-tracked  to  time  tn,  where  evaluation  of  w, l  is 
straight-forward.  Algorithmically,  this  is  implemented  by  evaluating  R  and  U  at  a  quadra¬ 
ture  point  xp  at  time  i„-i,  then  tracking  the  point  xp  from  tn_1  to  xp  =  X(<„ ;  xp,  t„_i)  at 
ip  and  determining  which  test  functions  are  nonzero  at  xp  at  tn  so  that  the  amount  of  mass 
associated  with  xp  can  be  added  to  the  corresponding  position  in  the  right-hand  side  vector 
in  the  global  discrete  linear  algebraic  system.  Notice  that  this  forward  tracking  has  no  effect 
on  the  solution  grid  or  the  data  structure  of  the  discrete  system.  Therefore,  the  forward 
tracking  algorithm  used  here  does  not  suffer  from  the  complication  of  distorted  grids,  which 
complicates  many  forward  tracking  algorithms,  and  is  a  major  attraction  of  the  backtracking 
in  characteristic  methods. 

4.2  Inflow  Boundaries 

If  fljj  <£_  Q,  either  E*,(/  intersects  the  inflow  boundary  rj/*,  or  T,nij  intersects  the  outflow 
boundary  F;,0) .  First,  consider  the  former  case  given  by  Equations  (3.5)  or  (3.6). 

The  first  term  on  the  left-hand  side  of  Equation  (3.5)  is  standard  in  finite  element  methods. 
The  second  terms  on  both  sides  are  standard  except  that  the  time  step  A £(^(x)  defined 
below  Equation  (2.6)  depends  on  x.  In  the  numerical  implementation,  we  calculate  these 
integrals  with  quadrature  points  at  time  tn.  Hence,  we  evaluate  A/ !/)  (x)  by  backtracking 
at  these  points.  For  each  quadrature  point  xp  E  Q ,:I  at  time  ,  we  need  to  track  the 
characteristic  X($;Xp,t„)  for  6  E  J„  to  determine  if  it  reaches  the  boundary  T,,  or  not.  If 
so,  we  calculate  the  time  G(xp)  when  the  characteristic  reaches  the  boundary  and  assign 
At/ (xp)  =  tn  —  ff  (xp);  otherwise,  At{l)  (xp)  =  At.  Notice  that  the  backtracking  algorithm 
is  used  only  to  calculate  A t^(x),  which  appears  in  the  diffusion-dispersion  term,  and  does 
not  affect  mass  conservation.  The  first  term  on  the  right-hand  side  of  Equation  (3.5)  can 
still  be  evaluated  by  a  forward  tracking  algorithm  as  in  Section  4.1. 

Notice  that  in  the  last  term  on  the  right-hand  side  of  Equation  (3.5)  ^^(x,  t)  is  defined 
at  the  space-time  boundary  I’-/' .  but  the  test  function  m;j  (x,  t)  =  w.;j(x,  tn),  where  x  = 
X(i„;  x,  t )  is  the  point  at  the  head  at  time  tn  corresponds  to  the  point  x  at  the  foot  at  time 
t.  Therefore,  we  use  a  forward  tracking  algorithm  to  calculate  this  term.  This  would  enforce 
the  integration  quadrature  at  the  space-time  boundary  F*/1  with  respect  to  a  fixed  spatial 


14 


grid  on  which  g':[]  (x,  t)  is  defined  and  track  forward  the  discrete  quadrature  points  chosen 
on  the  fixed  grid  at  the  space-time  boundary  Tjh  in  a  regular  fashion  to  time  tn,  where  one 
evaluates  wtj. 


Except  for  the  last  term  on  its  left-hand  side,  the  terms  in  Equation  (3.6)  are  similar  to 
those  in  Equation  (3.5).  The  evaluation  of  this  term  is  the  same  as  that  for  the  last  term 
on  the  right-hand  side  of  Equation  (3.5)  except  that  one  needs  to  use  forward  tracking  to 
evaluate  both  V U  and  the  test  function  Wjj. 


4.3  Outflow  Boundaries 

Consider  the  case  when  Q ([_  fi  and  E „  ,:I  intersects  the  outflow  boundary  E 1 .  The  nu¬ 
merical  scheme  is  given  by  Eciuations  (3.8)-(3.10).  We  only  discuss  the  evaluation  of  the 
integrals  in  Equations  (3.8)  and  (3.9)  since  the  evaluation  of  the  integrals  in  Equation  (3.10) 
is  similar. 

The  first  term  on  the  left-hand  side  of  Equation  (3.8)  is  standard.  The  second  terms  on 
both  sides  of  Equation  (3.8)  can  be  calculated  as  in  Section  4.2.  Keep  in  mind  that  the 
integrals  are  local  even  though  they  are  expressed  as  the  integrals  on  fi  at  time  Hence, 
At{ (x)  =  At  except  at  the  corner  of  =  (a,  b)  x  (c,  d)  where  T 1  and  T;/1  intersect.  The 
fact  that  inflow  and  outflow  boundaries  can  intersect  in  multiple  spatial  dimensions  makes 
the  implementation  more  complicated  than  that  for  one-dimensional  problems  where  inflow 
and  outflow  boundaries  do  not  meet  (as  long  as  the  one-dimensional  velocity  field  keeps  a 
definite  sign).  As  a  result  in  evaluating  the  second  terms  on  both  sides  of  Equation  (3.8),  we 
need  to  use  a  backward  tracking  algorithm  to  calculate  At{ (x)  near  the  corner  of  f 1  where 
the  inflow  boundary  Tj/*  and  outflow  boundary  V\(I>]  meet. 

In  Equation  (3.8),  the  four  integrals  defined  on  V 1  (with  the  first  V 1  integral  given 
by  Equation  (3.9))  are  standard  since  both  the  trial  function  U  and  the  test  functions 
Wjj  are  defined  on  EjW.  We  would  enforce  the  integration  quadrature  on  T<0).  Recall 
that  the  factor  Afl°*(x,  t)  in  some  of  these  terms  are  defined  by  (below  Equation  (2.6)) 
Afl0)(x,  t)  =  t  —  tn_i  except  when  the  characteristic  X($;x,  t)  meets  T;/1 .  In  this  case 
At(0)(x,  t)  =  t  —  t*(x,t)  where  t*(x,t)  E  Jn  is  the  time  when  X($;x,  t)  intersects  T;/1 .  In  the 
numerical  implementation,  we  simply  let  A fl0)(x,t)  =  t  —  tn_1,  except  near  the  corner  where 
the  inflow  boundary  and  the  outflow  boundary  T ^  intersect.  At  the  corner  region, 
we  use  a  backward  tracking  algorithm  to  locate  t*(x, t)  and  let  Afl0)(x, t)  =  t  —  t*(x, t). 
As  mentioned  in  Section  4.2,  the  use  of  backtracking  in  the  calculation  of  At* (x)  and 
Afl0)(x,  t)  does  not  effect  mass  conservation. 

The  first  term  on  the  right-hand  side  of  Equation  (3.8)  can  be  evaluated  by  a  forward 
tracking  algorithm  as  in  Sections  4. 1-4.2.  However,  notice  that  at  each  quadrature  point 
xp  E  f 2;.j  at  time  tn- i,  the  characteristic  X($;  xJ)?  t„_i)  may  intersect  Tj,0*.  In  the  current 
context,  we  need  to  use  a  forward  tracking  to  determine  if  X(0;  xp,  t„_i)  will  or  will  not 
intersect  Tj,0*.  In  the  latter  case  we  evaluate  Wij(xp,t,n)  as  in  Sections  4. 1-4.2.  In  the  former 
case,  we  need  to  locate  the  head  of  the  characteristic  at  the  space-time  boundary  T 1  and 
calculate  the  values  of  wl}  at  T 1  on  which  they  are  defined. 
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5  Description  of  Some  Other  Numerical  Methods 

In  this  section  we  briefly  describe  the  Galerkin  finite  element  method  (GAL),  the  quadratic 
Petrov-Galerkin  method  (Christie  et  al.  1976,  Barrett  and  Morton  1984,  and  Celia  et  al. 
1989),  the  cubic  Petrov-Galerkin  method  (Westerink  and  Shea  1989,  Bouloutas  and  Celia 
1991),  and  the  streamline  diffusion  finite  element  method  (Hughes  and  Brooks  1979,  Brooks 
and  Hughes  1982,  Hughes  and  Mallet  1986,  Hughes  1995,  Hansbo  and  Szepessy  1990,  Hansbo 
1992,  Johnson  and  Navert  1981,  Johnson  and  Pitkaranta  1986,  Johnson  1987,  Johnson  et  al. 
1990,  Johnson  and  Szepessy  1995,  Zhou  92,  and  Zhou  95).  For  simplicity  of  representation, 
we  assume  that  in  Equation  (2.1)  i?(x,  t)  =  1.  The  GAL,  QPG,  and  CPG  schemes  can  be 
unified  as  follows: 


Un  Wij  dx  +  XAt  |  Vwij  ■  (D„Vt/„)  dx  -  V„Un  ■  VwtJ  dxj 
=  Un- 1  Wjj  dx  -  (1  -  A) At  |  Vm;j  •  (D„_i  Vt/„_i)  dx 

-  jn  VnUn  ■  Vwij  dxj  +  At  |a  fn  wi:j  dx  +  (1  -  A)  /„_ i  WtJ  dx 
+  boundary  terms. 


(5.1) 


Here  A  G  [0, 1]  is  the  weighting  parameter  between  the  time  levels  t„_i  and  t„.  In  particular, 
A  =  1  and  0.5  yield  the  backward-Euler  (BE-)  and  the  Crank-Nicholson  (CN-)  schemes,  re¬ 
spectively.  The  trial  function  space  consists  of  the  standard  continuous  and  piecewise-bilinear 
polynomials.  The  test  functions  are  also  in  the  tensor  product  form  W;j(x,y)  =  Wj(x)w j(y). 
In  the  GAL  scheme,  W;(x)  and  Wj(y)  are  the  standard  one-dimensional  hat  functions.  In  the 
QPG  w,(x)  and  Wj(y)  are  constructed  by  adding  an  asymmetric  perturbation  to  the  original 
piecewise-linear  hat  functions 


w.i(x)  = 


X  -  Xi_  1  (x  -  xi_1)(xi  -  x) 

~^r  +  v — <A^ — ■ 

Xi+1  —  X  (x  —  Xj)(Xj+i  —  x) 


Ax 


—  v- 


{Axy 


0, 


,  X  G  [x,; ,  m 7; -j- 1  ] , 

otherwise  . 


(5.2) 


Here  v  =  3[coth(^|p)  —  for  constant  V  and  D.  For  variable  V  and  U,  one  replaces  A  by 
its  mean  on  each  element.  A  typical  one-dimensional  QPG  test  function  is  sketched  in  Figure 
1  (b).  As  defined  above,  the  two-dimensional  QPG  test  function  wrj(x,y)  is  just  a  tensor 
product  of  the  two  one-dimensional  QPG  test  functions  w,  (x)  and  W  j(y).  The  CPG  method 
was  derived  for  the  Crank-Nicholson  time  discretization.  In  the  CPG  w,(x)  and  w.t  (y)  are 
defined  as  the  original  piecewise  linear  hat  functions  with  a  symmetric  cubic  perturbation 
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added  to  each  nonzero  piece 


Wi(x)  =  l 


X  X;  I  (x  —  Xi_i)(Xi  —  x)(Xi_i  +  —  2x) 

+  7 - (Ax)3 - ’  A’  I  • 

(x  -  X;)(xt+1  ~  x)(X;  +  X;  .  ,  -  2x) 

7  /a  „vq  5  X  ^  mo  Pi+lJ) 


Ax 

X{-\- 1  x 


Ax 


(Ax)' 


(5.3) 


0. 


otherwise  . 


Here  7  =  5(7 u2  with  Cu  =  being  the  Courant  number.  For  variable  V  one  replaces  V 
by  its  arithmetic  mean  on  each  element.  A  typical  one-dimensional  CPG  test  function  is 
plotted  in  Figure  1  (c). 


The  SDM  is  a  type  of  discontinuous  Galerkin  FEM  and  applies  to  a  nonconservative 
analogue  of  Equation  (2.1).  For  the  following  nonconservative  advection-dispersion  equation 


£qu(x,  t)  =  ut  +  V(x,  t)  ■  Vtt  —  Y(D(x,  t)Vu)  = /(x,  t),  (x,  /)  G  E, 

u(x,  0)  =  tt0(x),  x  G  Q,  (2-1') 

u(x,  t)  =  0,  (x,  t)  G  r, 


the  discontinuous  trilinear  SDM  reads  as  follows:  find  a  piecewise-trilinear  (linear  in  time) 
function  U(x,t )  on  the  space-time  slab  E„  =  Q  X  ./„ .  which  is  discontinuous  in  time  at  1 
and  tn  and  satisfies  the  homogeneous  Dirichlet  boundary  condition,  such  that 


L  [Ut+V-  VC/]  [W  +  S(Wt  +  V  •  VIE)]  dxdt  +  L  VIE  •  (DVC/)  dxdt 

J  2jn  J 

-6  [  V  •  (DVU)(Wt  +  V  •  VIE)  dxdt  +  (  CZ+.^+.i  dx  (5.4) 

«/  /  -1  «/  u  u 

=  /s  /[ W  +  S(W,  +  V  .  VH/)]  dxdt  4-  Ja  U~_,  W,U  dx, 


for  any  test  function  IE  with  the  same  form  as  U.  Here  W^_1  =  lim  w(x,t)  and  Wn_1  = 


t — yt' 


lim  w(x,  t),  UQ  =  u o(x),  and  S  is  typically  chosen  to  be  of  0(h)  with  h  being  the  diameter 


of  the  space-time  partition  on  the  slab  E„.  The  third  term  on  the  left-hand  side  is  carried 
out  element-wise,  since  it  is  not  well-defined  for  piecewise-trilinear  functions. 


The  choice  of  S  has  significant  effects  on  the  accuracy  of  the  numerical  solutions.  If  <5 
is  chosen  too  small,  the  numerical  solutions  will  exhibit  oscillations.  If  <5  is  too  big,  the 
SDM  will  seriously  damp  the  numerical  solutions.  Unfortunately,  an  optimal  choice  of  S  is 
not  clear  and  is  heavily  problem-dependent.  Extensive  research  has  been  conducted  on  the 
SDM,  including  proper  choices  of  S  (Hughes  and  Brooks  1979,  Brooks  and  Hughes  1982, 
Hughes  and  Mallet  1986,  Hughes  1995,  Hansbo  and  Szepessy  1990,  Hansbo  1992,  Johnson 
and  Navert  1981,  Johnson  and  Pitkaranta  1986,  Johnson  1987,  Johnson  et  al.  1990,  Johnson 
and  Szepessy  1995,  Zhou  92,  and  Zhou  95).  Since  the  theme  of  this  paper  is  not  on  the 
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development  of  the  SDM,  in  our  numerical  experiments  we  use  a  generally  accepted  choice 
of  S  which  may  not  be  best  possible  for  a  given  problem.  According  to  Johnson  1987,  Johnson 
et  al.  1990,  and  Zhou  92,  we  set 


if  the  mesh  Peclet  number  |V| h  >  |D|  and  S  =  0  otherwise.  K  is  typically  to  be  1  or  0.5. 
In  our  numerical  experiments,  we  will  use  these  values  along  with  several  others  to  indicate 
the  general  behavior.  Moreover,  the  SDM  generally  increases  the  dimension  of  the  problem 
by  one  (although  the  measure  in  this  dimension  is  small).  For  Problem  (2.1’)  which  is  two- 
dimensions  in  space,  Equations  (5.4)  are  defined  on  three-dimensional  space-time  domain 
E„ .  Numerically,  one  has  to  partition  the  three-dimensional  “thick  slices”  into  tetrahedra  or 
prisms.  Usually  this  will  double  the  number  of  unknowns  in  GAL,  QPG,  CPG,  and  ELLAM 
schemes. 

While  the  SDM  can  capture  a  jump  discontinuity  of  the  exact  solution  in  a  thin  region, 
the  numerical  solution  may  develop  over-  and  under-slioots  about  the  exact  solution  within 
this  layer.  A  modified  SDM  with  improved  shock-capturing  properties  was  proposed  in 
[45,  49,  50],  which  consists  of  adding  a  “shock-capturing”  term  to  the  diffusion  by  introducing 
a  “cross-wind”  control  that  is  close  to  the  steep  fronts  or  “shocks”.  This  modified  SDM 
scheme  performs  much  better  in  terms  of  catching  the  steep  fronts  or  the  jump  discontinuities 
of  the  exact  solutions;  however,  it  leads  to  a  nonlinear  scheme  even  though  the  underlying 
governing  PDE  is  linear  and  involves  another  undetermined  parameter.  Thus,  we  will  not  use 
this  scheme  in  our  comparison  and  just  remind  the  reader  that  in  particular  cases  the  SDM 
may  perform  better  than  those  shown  in  the  examples  here  if  the  appropriate  modifications 
and  optimization  schemes  are  used. 

6  Computational  Results 

In  this  section  we  present  one-  and  two-dimensional  numerical  experiments  to  investigate 
the  performance  of  the  ELLAM  scheme  developed  in  this  paper  and  to  compare  it  with 
the  numerical  methods  described  in  Section  5.  The  numerical  experiments  contain  both 
examples  (with  analytical  solutions)  that  are  either  smooth  or  have  steep  fronts. 

6.1  The  One-Dimensional  Transport  of  a  Diffused  Square  Wave 

To  observe  the  performance  of  all  the  methods  in  Section  5  and  ELLAM  scheme  for  problems 
with  analytical  solutions  that  have  a  steep  front,  this  example  considers  the  transport  of  a 
one-dimensional  diffused  square  wave.  The  initial  condition  uq(x)  is  given  by 

1,  if C  xr\  C  (a,  6), 

u0(x)  =  (6.1) 

0,  otherwise. 

We  assume  that  the  one-dimensional  transport  equation  has  constant  coefficients  so  that  we 
can  find  the  analytical  solution  in  a  closed  form.  Homogeneous  inflow  and  outflow  Dirichlet 
boundary  conditions  are  specified  at  x  =  a  and  x  =  b.  As  long  as  the  diffused  square 
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wave  does  not  intersect  the  outflow  boundary  during  the  time  interval  [0,  T],  the  analytical 
solution  a (4,1)  can  be  expressed  as 


l  f — 

u(x,  t)  =  .  /  un(x  —  Vi  —  s )  exp  -  ds 

V  ’  ’  VImJ-oc  y  P  \  ADt  J 


erf 


'  x  —  Vi  —  xf 


~  erf 


'  x  —  Vi  —  xr 


(6.2) 


•JWt  J  V  Vwt 

where  erf(x)  =  ff  exp(— s2)ds  is  the  error  function. 


In  the  numerical  experiments  the  data  are  chosen  as  follows:  The  space  domain  is  (a,  b)  = 
(0,2),  the  time  interval  [0,T]  =  [0,1],  R  =  1,  V  =  1,  D  =  10-4,  and  /  =  0.  In  Equations 
(6.1)  and  (6.2),  Xi  =  0.2  and  xr  =  0.7,  so  that  the  diffused  square  wave  essentially  vanishes 
at  the  outflow  boundary  x  =  b  during  the  time  period  of  [0,T],  The  grid  size  Ax  =  pA 
is  chosen  so  that  the  analytical  solution  can  be  represented  properly.  The  backward-Euler 
Galerkin  linear  finite  element  (BE-GAL),  quadratic  Petrov-Galerkin  finite  element  (BE- 
QPG),  and  cubic  Petrov-Galerkin  finite  element  (BE-CPG)  solutions  are  plotted  against 
the  analytical  solution  in  Figures  2  (a)-(c)  for  At  =  Ajg,  gA,  and  Aqq,  respectively.  The 
ELL  AM  solution  is  also  plotted  in  Figure  2  (a)  for  At  =  A,  which  gives  a  Courant  number 
10  and  a  Peclet  number  100.  The  Crank-Nicholson  Galerkin  linear  finite  element  (CN- 
GAL),  quadratic  Petrov-Galerkin  finite  element  (CN-QPG),  and  cubic  Petrov-Galerkin  finite 
element  (CN-CPG)  solutions  are  plotted  in  Figures  3  (a)-(c)  for  At  =  jAj,  and 
respectively.  To  view  the  numerical  solutions  clearly,  we  do  not  plot  the  analytical  solution 
in  these  figures,  One  can  compare  the  CN-GAL,  CN-QPG,  and  CN-CPG  numerical  solutions 
with  the  analytical  one  in  Figures  2  (a)-(c).  The  streamline-diffusion  finite  element  (SDM) 
solutions  are  plotted  in  Figures  4  (a)-(b)  for  At  =  and  At  =  -Ay  The  SDM  solution  is 
also  plotted  in  Figure  4  (c)  for  Ax  =  At  =  A  to  further  observe  the  effect  of  the  choice  of 
the  parameter  <5. 


It  is  observed  that  the  BE-GAL,  BE-QPG,  and  BE-CPG  schemes  generate  almost  identi¬ 
cal  numerical  solutions.  With  a  time  step  of  At  =  ^A_ ,  |})e  backward  Euler  schemes  generate 
over-damped  numerical  solutions  without  any  observable  overshoot  or  undershoot.  As  the 
time  step  At  is  reduced  to  gA.  anq  _L_?  the  numerical  dispersion  is  reduced  considerably 
and  the  numerical  solutions  are  quite  close  to  the  analytical  solution.  With  a  time  step 
of  At  =  jAj5  the  CN-GAL  and  CN-QPG  solutions  display  overshoot  and  undershoot.  The 
maximum  and  minimum  values  of  CN-GAL  and  CN-QPG  solutions  are  1.212,  1.153,  —0.219, 
and  —0.153,  respectively.  The  CN-CPG  solution  also  has  many  wiggles  but  with  a  much 
smaller  magnitude  (Its  maximum  and  minimum  values  are  1.035  and  —0.031).  As  the  time 
step  At  is  reduced  to  ^Aj,  the  undershoot  and  overshoot  of  the  CN-GAL  solution  are  reduced 
by  about  40%  (The  maximum  and  minimum  values  are  1.132  and  —0.134).  The  undershoot 
and  overshoot  of  CN-QPG  solution  are  reduced  by  70%  (the  maximum  and  minimum  values 
are  1.047  and  —0.047)  and  are  comparable  to  those  of  CN-CPG  solution  (whose  maximum 
and  minimum  values  are  1.033  and  —0.033).  As  the  time  step  At  decreases  to  jA^,  the 
undershoot  and  overshoot  of  CN-GAL  solution  are  further  reduced  but  those  of  CN-QPG 
and  CN-CPG  solutions  do  not  change  much.  In  essence,  when  the  time  step  is  relatively 
large  (the  Courant  number  is  up  to  one),  the  CN-CPG  scheme  yields  better  solutions  than 
the  CN-GAL  and  CN-QPG  schemes. 
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With  a  time  step  At  =  the  SDM  solution  starts  to  approximate  the  analytical  solution. 
As  the  K  in  <5  decreases  from  1  to  0.001,  the  smearing  in  the  numerical  solutions  is  reduced 
considerably  and  the  overshoot/undershoot  is  also  reduced  slightly  (from  1.0952  and  —0.0577 
to  1.0714  and  —0.0721).  Thus,  the  optimal  value  of  K  seems  to  be  0.001  (the  smallest  of 
the  three  K  values).  As  the  time  step  At  decreases  to  pD,  the  SDM  solutions  become 
more  accurate  and  have  much  less  damping.  But  in  this  case  some  wiggles  appear  near  the 
locations  where  the  analytical  solution  has  steep  fronts.  Reducing  K  in  5  from  1  to  0.001 
gives  a  smaller  L2  error  in  the  SDM  solution  but  increases  the  overshoot/undershoot  slightly 
(from  1.0493  and  —0.0493  to  1.0592  and  —0.0592).  In  Figure  4  (c),  the  magnitude  of  the 
overshoot  and  undershoot  in  SDM  solutions  is  almost  doubled  (from  1.0520  and  —0.0548 
to  1.0932  and  —0.0952)  when  K  is  reduced  from  1  to  0.001.  Unfortunately,  there  is  no  a 
universal  rule  on  the  choice  of  the  K.  As  we  mentioned  above,  the  modified  SDM  with  a 
“shock  capturing”  property  should  generate  better  numerical  solutions  than  those  shown 
here.  However,  one  has  to  solve  a  nonlinear  system  even  though  the  underlying  PDE  is  a 
linear  one,  and  must  face  choosing  an  additional  undetermined  parameter.  In  contrast,  with 
a  fairly  large  time  step  At  =  j-}  the  ELLAM  scheme  yields  a  very  accurate  numerical  solution 
that  is  better  than  any  one  of  the  GAL,  QPG,  CPG,  and  SDM  solutions  with  even  much 
finer  time  steps.  Moreover,  the  ELLAM  scheme  uses  only  half  the  number  of  unknowns  as 
in  the  SDM  and  does  not  require  optimization  of  indefinite  constants. 


6.2  A  Two-Dimensional  Rotating  Gaussian  Pulse 

This  example  considers  the  transport  of  a  two-dimensional  rotating  Gaussian  pulse.  The 
spatial  domain  is  Q,  =  (—0.5,  0.5)  x  (— 0.5, 0.5),  the  rotating  field  is  imposed  as  V\(x,  y)  =  —4 y, 
and  V2 (x,y)  =  Ax.  The  time  interval  is  [0,T]  =  [0,7t/2],  which  is  the  time  period  required 
for  one  complete  rotation.  The  initial  condition  u0(x,y )  is  given  by 

^)=exp  (6.3) 


where  xc,  yc,  and  a  are  the  centered  and  standard  deviations,  respectively.  The  corresponding 
analytical  solution  for  Equation  (2.1)  with  R  =  1,  a  constant  diffusion  coefficient  D,  and 
/  =  0  is  given  by 


"(4//-0  = 


2  a2 


exp 


(x  -  xcf  +  (y  -  yc): 


2  a2  +  ADt  y  2cr2  +  ADt 

where  x  =  .tcos(4<)  +  ysin(4i)  and  y  =  — ®sin(4f)  +  ycos(At). 


(6.4) 


In  the  numerical  experiments,  the  data  are  chosen  as  follows:  D  =  10-4,  xc  =  —0.25, 
yc  =  0,  a  =  0.0447  which  gives  2cr2  =  0.0040.  A  uniform  spatial  grid  of  Ax  =  Ay  =  -p  is 
used  (in  all  the  plots  and  most  of  the  experiments  in  Table  1),  unless  it  is  specified  otherwise, 
in  which  case  a  uniform  spatial  grid  of  Ax  =  Ay  =  ^  is  used.  It  is  easy  to  see  that  the 
u0(x,  y)  defined  by  (6.3)  is  centered  at  {xr.  y,.)  with  a  minimum  value  0  and  a  maximum 
value  1.  Its  surface  and  contour  plots  are  presented  in  Figures  5  (a)-(b).  Figures  5  (c)-(d) 
are  the  surface  and  contour  plots  for  the  analytical  solution,  which  has  a  minimum  value  0 
and  a  maximum  value  0.8642  (due  to  the  effect  of  diffusion). 
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This  problem  provides  an  example  for  a  homogeneous  two-dimensional  advection-dispersion 
equation  with  a  variable  velocity  field  and  a  known  analytical  solution.  Moreover,  this  prob¬ 
lem  changes  from  the  advection  dominance  in  most  of  the  domain  to  the  diffusion  dominance 
in  the  region  that  is  close  to  the  origin.  These  types  of  problems  often  arise  in  many  important 
applications  and  are  more  difficult  to  simulate  compared  with  purely  advection-dominated 
problems.  This  example  has  been  used  widely  to  test  for  numerical  artifacts  of  different 
schemes,  such  as  numerical  stability  and  numerical  dispersion,  spurious  oscillations,  and 
phase  errors. 

In  our  experiments,  we  have  systematically  varied  the  time  steps  to  examine  the  perfor¬ 
mance  of  each  method,  using  a  uniform  spatial  grid  of  Ax  =  Ay  =  ^  in  most  experiments. 
This  is  because  the  temporal  errors  dominate  the  numerical  solutions  with  all  the  methods 
other  than  the  ELL  AM  schemes.  In  this  case  the  maximum  grid  Peclet  number  reaches  442. 
The  grid  Peclet  number  at  the  center  of  the  Gaussian  Pulse  is  about  156.  We  also  perform 
some  experiments  with  a  finer  spatial  grid  of  Ax  =  Ay  =  In  this  case  the  maximum  grid 
Peclet  number  reaches  295  and  the  grid  Peclet  number  at  the  center  of  the  Gaussian  Pulse 
is  about  104.  All  comparative  methods  tested  yield  strongly  non-symmetric  systems  while 
the  ELLAM  scheme  inherently  symmetrizes  its  discrete  algebraic  system.  We  use  a  precon¬ 
ditioned  conjugate  gradient  square  algorithm  (PCGS)  to  solve  these  systems  even  though 
this  places  ELLAM  at  a  disadvantage.  In  Table  1  we  present  the  minimum  and  maximum 
values  of  the  numerical  solutions  with  each  method,  and  the  per  time  step  CPU  and  the 
overall  CPU  each  method  consumed,  which  was  measured  on  a  SGI  Indigo  Workstation.  We 
realize,  of  course,  that  some  code  optimization  may  be  possible  but  feel  that  these  timings 
are  representative  of  each  scheme’s  efficiency  on  these  model  problems.  The  surface  and 
contour  plots  for  selected  runs  of  each  method  in  Table  1  are  presented  in  Figures  5  -  12. 

6.2.1  The  ELLAM  Simulation 

The  ELLAM-Euler  solution  is  obtained  by  using  a  time  step  of  At  =  in  solving  the  ELLAM 
scheme  and  using  an  Euler  quadrature  with  a  micro-time  step  of  A tf  =  in  tracking  the 
characteristics.  The  maximum  Courant  number  reaches  14,  while  the  Courant  number  at 
the  center  of  the  Gaussian  Pulse  is  5.  The  ELLAM-Euler  solution  has  a  minimum  value  0 
and  a  maximum  value  0.8302,  whose  surface  and  contour  plots  are  given  in  Figures  6  (a)-(b). 
As  one  can  see  from  these  plots  and  Table  1,  a  very  accurate  numerical  solution  is  obtained 
in  about  5  and  a  half  minutes.  The  ELLAM-RK  solution  is  obtained  by  using  a  time  step  of 
At  =  -^  in  solving  the  ELLAM  scheme  and  using  a  second-order  Runge-Kutta  quadrature 
with  a  micro-time  step  of  A  tf  =  ^-  in  tracking  the  characteristics.  In  this  case,  the  maximum 
Courant  number  reaches  19,  while  the  Courant  number  at  the  center  of  the  Gaussian  Pulse  is 
about  6.7.  The  ELLAM-RK  solution  has  a  minimum  value  0  and  a  maximum  value  0.8630, 
whose  plots  are  in  Figures  6  (c)-(d).  The  use  of  a  more  accurate  second-order  Runge-Kutta 
tracking  algorithm  enables  us  to  significantly  reduce  the  number  of  micro-time  steps  (from 
80  in  an  Euler  tracking  to  4  in  the  Runge-Kutta  tracking)  in  tracking  characteristics  and 
so  the  CPU  time  per  global  time  step  (from  17  seconds  in  ELLAM-Euler  to  10.2  seconds  in 
ELLAM-RK).  Moreover,  the  number  of  global  time  steps  is  reduced  from  20  in  the  ELLAM- 
Euler  simulation  to  15  in  the  ELLAM-RK  simulation.  Thus,  the  ELLAM-RK  simulation 
further  reduces  the  overall  CPU  time  to  2  minutes  and  33  seconds. 
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6.2.2  The  BE-GAL  and  BE-QPG  Simulation 


Due  to  its  unconditional  stability  and  simplicity  in  implementations,  the  fully  implicit  back¬ 
ward  Euler  temporal  discretization  still  dominates  most  production  codes  in  engineering 
applications.  Thus,  we  carry  out  extensive  experiments  to  investigate  the  performance  of 
this  discretization.  With  a  time  step  of  At  =  which  gives  a  maximum  Courant  number 
of  2.84  and  a  Courant  number  of  1  at  the  center  of  the  Gaussian  Pulse,  the  BE-GAL  and 
BE-QPG  solutions  have  minimum  values  of  0  and  —0.0002  and  maximum  values  of  0.2546 
and  0.2003,  respectively.  Recall  that  the  exact  solution  has  a  maximum  value  of  0.8642,  the 
BE-GAL  and  BE-QPG  solutions  are  excessively  over-damped.  Moreover,  the  BE-GAL  and 
BE-QPG  schemes  require  more  iterations  in  the  PCGS  solver  than  the  ELLAM  does,  because 
they  yield  strongly  non-symmetric  coefficient  matrices.  The  BE-GAL  and  BE-QPG  solutions 
with  a  much  finer  time  step  of  At  =  ^  are  presented  in  Figures  7  (a)-(d).  The  minimum 
values  are  0  for  BE-GAL  solution  and  —0.0003  for  BE-QPG  solution,  while  the  maximum 
values  are  0.4517  for  BE-GAL  solution  and  0.3486  for  BE-QPG  solution.  These  solutions 
are  still  very  diffusive  and  are  considerably  deformed,  especially  for  the  BE-QPG  solution. 
The  more  severe  deformation  in  the  BE-QPG  solution  is  due  to  the  effect  of  grid  orientation 
incurred  by  the  upwinding  in  the  QPG  method  (refer  to  Ewing  1984).  With  the  same  time 
step  of  At  =  we  also  reduce  the  spatial  grid  from  Ax  =  Ay  =  A  to  Ax  =  Ay  =  A 
to  observe  the  improvement  of  the  numerical  solutions.  The  BE-QPG  solution  has  a  slight 
improvement  while  the  BE-GAL  has  essentially  no  improvement.  However,  the  CPU  time 
has  been  significantly  increased.  With  a  comparable  overall  CPU  time  we  could  use  a  much 
finer  time  step  of  At  =  and  still  use  the  coarse  spatial  grid  of  Ax  =  Ay  =  A.  in  this 
case  the  numerical  solutions  have  more  visible  improvement.  This  shows  that  even  with  a 
time  step  of  At  =  ^  and  a  spatial  grid  of  Ax  =  Ay  =  A?  the  temporal  dominance  still 
dominates  the  numerical  solutions  in  the  backward  temporal  discretization.  Note  that  for 
Ax  =  Ay  =  A  and  A/  =  the  corresponding  maximum  Courant  number  and  the  Courant 
number  at  the  center  of  the  Gaussian  Pulse  are  0.7  and  0.25,  respectively. 

To  obtain  BE-GAL  and  BE-QPG  solutions  with  reasonable  accuracy,  we  proceed  further. 
With  a  time  step  of  At  =  ^q,  we  reduce  the  spatial  grid  from  Ax  =  Ay  =  A  to  Ax  = 
Ay  =  A.  Again  one  sees  slight  improvement  in  the  BE-QPG  solution  and  no  improvement 
in  the  BE-GAL  solution,  but  a  significant  increase  in  the  overall  CPU  time.  Using  less 
overall  CPU  time,  we  could  use  the  original  coarse  spatial  grid  of  Ax  =  Ay  =  A  but  a  finer 
time  step  of  At  =  Note  that  for  Ax  =  Ay  =  A  and  At  =  %  the  corresponding 

maximum  Courant  number  and  the  Courant  number  at  the  center  of  the  Gaussian  Pulse  are 
0.28  and  0.1,  respectively.  This  shows  that  the  temporal  error  still  dominates  the  BE-GAL 
and  BE-QPG  solutions.  Our  last  numerical  experiments  with  the  backward  Euler  temporal 
discretization  used  a  spatial  grid  of  Ax  =  Ay  =  A  and  a  time  step  of  At  =  g^.  The 
minimum  values  are  0  for  BE-GAL  solution  and  —0.0017  for  BE-QPG  solution,  while  the 
maximum  values  are  0.7401  for  BE-GAL  solution  and  0.5554  for  BE-QPG  solution.  Their 
surface  and  contour  plots  are  presented  in  Figures  8  (a)-(d).  Note  that  in  this  case  the 
maximum  Courant  number  and  the  Courant  number  at  the  center  of  the  Gaussian  Pulse  are 
0.09  and  0.03,  respectively.  However,  the  BE-GAL  solution  is  still  not  comparable  with  the 
two  ELLAM  solutions.  The  BE-QPG  solution  is  even  much  worse.  In  fact,  Figure  8  (c)-(d) 
show  that  the  BE-QPG  solution  has  severe  deformation.  However,  the  overall  CPU  time 
is  more  than  11  hours  for  the  BE-GAL  solution  and  more  than  15  hours  for  the  BE-QPG 
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solution.  This  is  in  contrast  to  the  4  minutes  and  39  seconds  for  the  ELLAM-Euler  solution 
and  2  minutes  33  seconds  for  the  ELLAM-RK  solution. 

Therefore,  even  though  the  backward  Euler  temporal  discretization  is  unconditionally 
stable  and  simple  to  implement,  extremely  small  time  steps  have  to  be  used  in  these  schemes, 
not  for  the  purpose  of  stability,  but  for  the  purpose  of  comparative  accuracy  (refer  to  Ewing 
1984).  Consequently,  this  significantly  reduces  the  efficiency  of  the  simulation. 

6.2.3  The  CN-GAL,  CN-QPG,  and  CN-CPG  Simulation 

The  CN-GAL  and  CN-QPG  solutions  are  presented  in  Figures  9  (a)-(d)  for  a  spatial  grid 
of  Ax  =  Ay  =  and  a  time  step  of  At  =  ^Jq.  The  CN-GAL  solution  has  minimum  and 
maximum  values  —0.1564  and  0.7861.  Severe  undershoot,  deformation,  and  phase  errors 
are  observed  in  the  CN-GAL  solution  in  Figures  9  (a)-(b).  The  CN-QPG  solution  has  a 
minimum  value  —0.0978  and  a  maximum  value  0.6197,  respectively.  Hence,  the  CN-QPG 
solution  has  about  40%  less  undershoot  than  the  CN-GAL  solution,  but  it  also  has  serious 
damping,  phase  error,  and  deformation.  The  CN-CPG  solution  is  not  available  (unbounded) 
for  the  time  step  and  spatial  grid.  Note  that  the  maximum  Courant  number  is  2.84  and 
the  Courant  number  is  1  at  the  center  of  the  Gaussian  Pulse  in  the  current  circumstances. 
Also,  note  that  the  Crank-Nicholson  temporal  discretization  yields  more  accurate  numerical 
solutions  than  the  backward  Euler  temporal  discretization,  due  to  its  higher-order  temporal 
accuracy.  The  overall  CPU  time  in  this  case  is  about  25  minutes  for  the  CN-GAL  solution 
and  about  30  minutes  for  the  CN-QPG  solution.  The  BE-GAL  and  BE-QPG  solutions  do 
not  have  undershoot,  but  the  CN-GAL  and  CN-QPG  solutions  do  exhibit  serious  problems  in 
this  regard  and  indicate  a  considerable  disadvantage  for  the  Crank-Nicholson  discretization. 

We  further  reduce  the  time  step  to  At  =  in  the  numerical  simulation.  In  this  case, 
the  CN-CPG  solution  is  also  available.  We  present  the  results  in  Table  1  and  Figures  10 
(a)-(d)  and  11  (a)-(b).  The  maximum  values  are  0.8438  for  the  CN-GAL  solution,  0.6412  for 
the  CN-QPG  solution,  and  0.8555  for  the  CN-CPG  solution,  while  the  minimum  values  are 
—0.0159  for  the  CN-GAL  solution,  —0.0081  for  the  CN-QPG  solution,  and  —0.0002  for  the 
CN-CPG  solution.  As  one  can  see,  the  numerical  solutions  have  been  improved  considerably 
and  the  undershoot  in  the  solutions  has  been  rapidly  reduced.  However,  these  solutions  still 
have  deformation,  especially  in  the  case  of  the  CN-QPG  solution. 

6.2.4  The  SDM  simulation 

The  surface  and  contour  plots  of  SDM  solutions  are  plotted  in  Figures  11  (c)-(d)  and  12 
(a)-(d)  for  a  time  step  of  At  =  and  Ax  =  Ay  =  -tj.  The  undetermined  parameter  K  in 
(5.5)  equals  to  0.5,  0.1,  and  0.001,  respectively.  As  K  decreases  from  0.5  to  0.1  and  then  to 
0.001,  the  maximum  and  minimum  values  of  the  corresponding  SDM  solutions  change  from 
0.7089  and  —0.0147  to  0.8250  and  —0.0021,  and  then  to  0.8281  and  —0.0019.  Namely,  the 
SDM  solutions  have  eliminated  almost  all  the  damping  and  undershoot  and  become  more 
accurate.  The  numerical  solutions  will  no  longer  improve  as  one  further  reduces  the  value 
of  K.  The  SDM  solutions  have  no  phase  error  or  deformation,  but  does  require  the  most 
CPU  time  per  time  step  since  it  has  double  the  number  of  unknowns  as  those  for  the  other 
methods.  This  in  turn  requires  more  iterations  in  solving  the  linear  system.  Furthermore, 
on  each  (space-time)  cell,  the  SDM  has  eight  basis  functions  which  are  the  tensor  product  of 
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three  univariate  functions,  while  all  other  methods  have  four  basis  functions  on  each  (space) 
cell  which  are  the  tensor  product  of  two  univariate  functions. 

In  summary,  one  sees  from  Table  1  and  Figures  5-12  that  the  ELLAM  scheme  is  the  most 
(CPU)  cost  effective  per  time  step  and  is  much  more  cost  effective  over  all,  since  ELLAM 
scheme  outperforms  the  other  methods  tested  with  much  fewer  time  steps. 

We  refer  readers  to  the  works  of  Ewing  and  Wang  (1991,  1993a,  1996),  Wang  (1992), 
Wang  and  Ewing  (1995),  and  Wang  et  al  (1995a)  for  the  asymptotic  convergence  rates  of 
the  ELLAM  schemes  in  space  and  time,  where  numerical  experiments  were  performed  and 
theoretical  convergence  error  estimates  were  derived  for  the  asymptotic  convergence  rates 
of  the  ELLAM  schemes  for  first-order  linear  hyperbolic  equations  and  advection-diffusion 
equations  in  one  space  dimension.  A  theoretical  error  estimate  for  an  ELLAM  scheme  for 
a  two-dimensional  advection-diffusion  equations  with  constant  coefficients  was  also  outlined 
in  Wang  (1992). 

7  Summary 

In  this  paper  we  have  developed  an  Eulerian-Lagrangian  localized  adjoint  method  for  two- 
dimensional  linear  advection-dispersion  equations  with  general  inflow  and  outflow  boundary 
conditions  based  on  the  approach  presented  in  Wang  (1992).  The  derived  numerical  scheme 
conserves  mass  and  treats  any  combinations  of  inflow  and  outflow  Diriclrlet,  Neumann,  and 
Robin  boundary  conditions  in  a  systematic  manner. 

Traditional  forward-tracked  characteristic  methods  or  particle  methods  advance  the  grids 
following  the  characteristics,  which  typically  result  in  severely  distorting  the  evolving  grids 
even  though  the  initial  grids  were  uniform.  This  greatly  complicates  the  solution  proce¬ 
dures.  Many  characteristic  methods  including  certain  ELLAM  schemes  have  been  developed 
using  a  backtracking  algorithm  to  avoid  these  problems  (Arbogast  et  al.  1992,  Arbogast 
and  Wheeler  1995,  Binning  1994,  Binning  and  Celia  1996,  Celia  et  al.  1990,  Dahle  et 
al.  1990,  Douglas  and  Russell  1982,  Espedal  and  Ewing  1987,  Ewing  1984,  Ewing  and 
Wang  1991, 1993a,  1993b,  1994, 1996,  Russell  1990,  Wang  et  al.  1995).  However,  for  multi¬ 
dimensional  problems,  backtracked  characteristic  methods  require  significant  effort,  due  to 
the  need  to  define  the  geometry  at  time  tn_  1  which  requires  the  tracking  of  points  along 
the  boundary  of  the  element  and  subsequent  interpolation  and  mapping  onto  the  fixed  spa¬ 
tial  grid  at  the  previous  time  level  i.  This  approach  is  computationally  very  intensive, 
especially  when  part  or  all  of  the  element  being  mapped  intersects  a  space-time  boundary 
(Binning  and  Celia  1996,  Celia  1994). 

The  ELLAM  scheme  in  this  paper  uses  a  forward-tracking  approach  (Dahle  et  al.  1995, 
Healy  and  Russell  1993,  Russell  and  Trujillo  1990,  Wang  1992)  to  track  quadrature  points 
at  time  tn_1  in  evaluating  the  storage  terms  and  inflow  boundary  terms  on  the  right-hand 
side  of  equations  (3.4),  (3.5),  (3.6),  and  (3.8).  Thus,  this  forward  tracking  scheme  has  no 
effect  on  the  underlying  grid  or  the  data  structure  of  the  discrete  system.  Furthermore,  the 
scheme  uses  a  backtracking  of  characteristics  to  evaluate  the  x-dependent  time  step  A (x) 
in  the  diffusion-dispersion  term. 
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In  this  paper  we  have  performed  one-  and  two-dimensional  numerical  experiments  to 
observe  the  performance  of  the  ELLAM  scheme  and  to  compare  it  with  many  intensely  in¬ 
vestigated  and  well  received  methods,  such  as  the  standard  Galerkin  finite  element  method, 
quadratic  Petrov-Galerkin  finite  element  method,  and  cubic  Petrov-Galerkin  method,  which 
use  either  a  backward  Euler  or  a  Crank-Nicholson  temporal  discretization,  as  well  as  the 
streamline  diffusion  finite  element  method.  The  numerical  experiments  show  that  the  EL¬ 
LAM  scheme  has  generated  very  accurate  numerical  solutions  (compared  with  the  other 
methods  considered)  even  though  a  much  larger  time  step  is  used  in  the  ELLAM  scheme. 
Consequently,  the  ELLAM  scheme  has  a  significantly  enhanced  efficiency.  In  the  context  of 
one-dimensional  first-order  linear  hyperbolic  equations,  a  more  extensive  comparison  of  the 
ELLAM  schemes  with  many  other  well-regarded  methods,  including  the  continuous  and  dis¬ 
continuous  Galerkin  finite  element  method  (Falk  and  Richter  1992,  Johnson  1987,  Johnson 
and  Pitkaranta  1986,  Richter  1988),  the  monotonic  upstream-centered  scheme  for  conserva¬ 
tion  laws  (Colella  1985,  van  Leer  1984)  and  the  essentially  non-oscillatory  scheme  (Dawson 
1991,  Harten  et  al.  1987,  Shu  and  Osher  1988),  as  well  as  the  methods  described  in  Section 
5,  can  be  found  in  the  work  of  Wang  et  al.  (Wang,  Al-Lawatia,  and  Telyakovskiy) . 

Finally,  we  point  out  that  the  Eulerian  methods  are  relatively  easy  to  formulate  and  to 
implement,  in  general.  In  contrast,  due  to  the  use  of  the  Lagrangian  coordinates,  character¬ 
istic  methods  (including  the  ELLAM  scheme)  typically  require  more  implementational  work, 
especially  for  multi-dimensional  problems. 

Acknowledgments: 

We  wish  to  acknowledge  the  many  helpful  discussions  and  suggestions  regarding  this  work 
which  were  provided  by  Professors  Michael  A.  Celia  and  Thomas  F.  Russell. 

References 

[1]  T.  Arbogast,  A.  Cliilakapati,  and  M.F.  Wheeler:  A  characteristic-mixed  method  for 
contaminant  transport  and  miscible  displacement,  Russell  et  al.  (eds.)  Computational 
Methods  in  Water  Resources  IX  1,  Computational  Mechanics  Publications  and  Elsevier 
Applied  Science,  London  and  New  York,  (1992),  77-84. 

[2]  T.  Arbogast  and  M.F.  Wheeler:  A  characteristic-mixed  finite  element  method  for 
advection-dominated  transport  problems,  SIAM  J.  Num.  Anal,  32,  (1995),  404-424. 

[3]  R.  Bank,  J.  Burger,  W.  Fichtner,  and  R.  Smith:  Some  upwinding  techniques  for  finite 
element  approximations  of  convection  diffusion  equations,  Num.  Math.,  58,  (1990),  185— 
202. 

[4]  A.M.  Baptista:  Solution  of  advection-dominated  transport  by  Eulerian-Lagrangian 
methods  using  the  backwards  method  of  characteristics.  Pli.D.  Thesis,  Massachusetts 
Institute  of  Technology,  1987. 

[5]  J.W.  Barrett  and  K.W.  Morton:  Approximate  symmetrization  and  Petrov-Galerkin 
methods  for  diffusion-convection  problems,  Comp.  Meth.  Appl.  Mech.  Engrg.,  45, 
(1984),  97-122. 


25 


[6]  J.P.  Benque  and  J.  Ronat:  Quelques  difficulties  des  modeles  numeriques  en  hydraulique, 
Comp.  Meth.  Appl.  Mech.  Engrg .,  Glowinski  and  Lions  (eds.),  North-Holland,  (1982), 
471-494. 

[7]  P.J.  Binning:  Modeling  unsaturated  zone  flow  and  contaminant  in  the  air  and  wa¬ 
ter  phases.  Pli.D.  Thesis,  Department  of  Civil  Engineering  and  Operational  Research, 
Princeton  University,  1994. 

[8]  P.J.  Binning  and  M.A.  Celia:  A  finite  volume  Eulerian-Lagrangian  localized  adjoint 
method  for  solution  of  the  contaminant  transport  equations  in  two-dimensional  multi¬ 
phase  flow  systems,  Water  Resource  Research,  32,  (1996),  103-114. 

[9]  E.T.  Bouloutas  and  M.A.  Celia:  An  improved  cubic  Petrov-Galerkin  method  for  simula¬ 
tion  of  transient  advection-diffusion  processes  in  rectangularly  decomposable  domains, 
Comp.  Meth.  Appl.  Mech.  Engrg.,  91,  (1991),  289-308. 

[10]  A.  Brooks  and  T.J.R.  Hughes:  Streamline  upwind  Petrov-Galerkin  formulations  for 
convection  dominated  flows  with  particular  emphasis  on  the  incompressible  Navier- 
Stokes  equations,  Comp.  Meth.  Appl.  Mech.  Engrg.,  32,  (1982),  199-259. 

[11]  M.A.  Celia,  I.  Herrera,  E.T.  Bouloutas,  and  J.S.  Kindred:  A  new  numerical  approach 
for  the  advective-diffusive  transport  equation,  Num.  Meth.  PDEs,  5,  (1989),  203-226. 

[12]  M.A.  Celia,  T.F.  Russell,  I.  Herrera,  and  R.E.  Ewing:  An  Eulerian-Lagrangian  localized 
adjoint  method  for  the  advection-diffusion  equation,  Advances  in  Water  Resources,  13, 
(1990),  187-206. 

[13]  M.A.  Celia  and  S.  Zisman:  An  Eulerian-Lagrangian  localized  adjoint  method  for  reactive 
transport  in  groundwater,  Gambolati  et  al.  (eds.),  Computational  Methods  in  Water 
Resources  VIII.  1,  Springer,  (1990),  383-392. 

[14]  M.A.  Celia  and  L.A.  Ferrand:  A  comparison  of  ELLAM  formulations  for  simulation  of 
reactive  transport  in  groundwater,  Advances  in  Hydro-Science  and  Engineering ,  1(B), 
Wang  (ed.),  University  of  Mississippi  Press,  (1993),  1829-1836. 

[15]  M.A.  Celia:  Eulerian-Lagrangian  localized  adjoint  methods  for  contaminant  transport 
simulations,  Peters  et  al.  (eds.),  Computational  Methods  in  Water  Resources  X,  1,  Water 
Science  and  Technology  Library,  12,  Kluwer  Academic  Publishers,  Dordrecht,  Nether¬ 
lands,  (1994),  207-216. 

[16]  I.  Christie,  D.F.  Griffiths,  A.R.  Mitchell,  and  O.C.  Zienkiewicz:  Finite  element  methods 
for  second  order  differential  equations  with  significant  first  derivatives,  Int.  J.  Num. 
Engrg.,  10,  (1976),  1389-1396. 

[17]  P.  Colella:  A  direct  Eulerian  MUSCL  scheme  for  gas  dynamics,  SIAM  J.  Sci.  Stat. 
Comp.  6,  (1985),  104-117. 

[18]  R.A.  Cox  and  T.  Nishikawa:  A  new  total  variation  diminishing  scheme  for  the  solution  of 
advective-dominant  solute  transport,  Water  Resource  Research,  27,  (1991),  2645-2654. 

[19]  H.K.  Dahle,  M.S.  Espedal,  R.E.  Ewing,  and  O.  Saevareid:  Characteristic  adaptive  sub- 
domain  methods  for  reservoir  flow  problems,  Num.  Meth.  PDEs ,  (1990),  279-309. 

[20]  H.K.  Dahle,  R.E.  Ewing,  and  T.F.  Russell:  Eulerian-Lagrangian  localized  adjoint  meth¬ 
ods  for  a  nonlinear  convection-diffusion  equation,  Comp.  Meth.  Appl.  Mech.  Engrg.,  122, 
(1995),  223-250. 


26 


[21]  C.N.  Dawson:  Godunov-mixed  methods  for  advective  flow  problems  in  one  space  di¬ 
mension,  SIAM  J.  Num.  Anal,  28,  (1991),  1282-1309. 

[22]  L.  Demkowitz  and  J.T.  Oden:  An  adaptive  characteristic  Petrov-Galerkin  finite  element 
method  for  convection-dominated  linear  and  nonlinear  parabolic  problems  in  two  space 
variables,  Comp.  Meth.  Appl.  Mech.  Engrg.,  55,  (1986),  63-87. 

[23]  J.  Douglas,  Jr.  and  T.F.  Russell:  Numerical  methods  for  convection-dominated  diffusion 
problems  based  on  combining  the  method  of  characteristics  with  finite  element  or  finite 
difference  procedures,  SIAM  J.  Num.  Anal.,  19,  (1982),  871-885. 

[24]  K.  Eriksson  and  C.  Johnson:  Adaptive  streamline  diffusion  finite  element  methods  for 
stationary  convection-diffusion  problems,  Math.  Comp.,  60,  (1993),  167-188 

[25]  M.S.  Espedal  and  R.E.  Ewing:  Characteristic  Petrov-Galerkin  subdomain  methods  for 
two-pliase  immiscible  flow,  Comp.  Meth.  Appl.  Mech.  Engrg.,  64,  (1987),  113-135. 

[26]  R.E.  Ewing,  T.F.  Russell,  and  M.F.  Wheeler:  Simulation  of  miscible  displacement  using 
mixed  methods  and  a  modified  method  of  characteristics,  SPE  12241,  (1983),  71-81. 

[27]  R.E.  Ewing  (ed.):  The  Mathematics  of  Reservoir  Simulation.  Frontiers  in  Applied  Math¬ 
ematics,  1,  SIAM,  Philadelphia,  1984. 

[28]  R.E.  Ewing:  Operator  splitting  and  Eulerian-Lagrangian  localized  adjoint  methods  for 
multiphase  flow,  The  Mathematics  of  Finite  Elements  and  Applications  VII,  Whiteman 
(ed.),  Academic  Press,  Inc.,  San  Diego,  CA,  (1991),  215-232. 

[29]  R.E.  Ewing  and  H.  Wang:  Eulerian-Lagrangian  localized  adjoint  methods  for  linear 
advection  equations,  Computational  Mechanics  ’91,  Springer  International,  (1991),  245- 
250. 

[30]  R.E.  Ewing  and  H.  Wang:  Eulerian-Lagrangian  localized  adjoint  methods  for  linear 
advection  or  advection-reaction  equations  and  their  convergence  analysis,  Comp.  Mech., 
12,  (1993a),  97-121. 

[31]  R.E.  Ewing  and  H.  Wang:  An  Eulerian-Lagrangian  localized  adjoint  method  for 
variable-coefficient  advection-reaction  problems,  Advances  in  Hydro- Science  and  En¬ 
gineering  1(B)  (Wang,  ed.),  University  of  Mississippi  Press,  (1993b),  2010-2015. 

[32]  R.E.  Ewing  and  H.  Wang:  An  Eulerian-Lagrangian  localized  adjoint  method 
with  exponential-along-cliaracteristic  test  functions  for  variable-coefficient  advective- 
diffusive-reactive  equations.  Choi  et  al.  (eds)  Proceedings  of  KAIST  Mathematical  Work¬ 
shop,  Analysis  and  Geometry,  8,  Taejon,  Korea,  (1993c),  77-91. 

[33]  R.E.  Ewing  and  H.  Wang:  Eulerian-Lagrangian  localized  adjoint  methods  for  variable- 
coefficient  advective-diffusive-reactive  equations  in  groundwater  contaminant  transport, 
Advances  in  Optimization  and  Numerical  Analysis,  Mathematics  and  Its  Applications, 
275,  Gomez  and  Hennart  (eds.),  Kluwer  Academic  Publishers,  Dordrecht,  Netherlands, 
(1994),  185-205. 

[34]  R.E.  Ewing  and  H.  Wang:  An  optimal-order  error  estimate  to  Eulerian-Lagrangian 
localized  adjoint  method  for  variable-coefficient  advection-reaction  problems.  SIAM  J. 
Num.  Anal,  33(1),  (1996),  318-348. 

[35]  R.S.  Falk  and  G.R.  Richter:  Local  error  estimates  for  a  finite  element  method  for  hy¬ 
perbolic  and  convection-diffusion  equations,  SIAM  J.  Num.  Anal,  29,  (1992),  730-754. 


27 


[36]  A.O.  Garder,  D.W.  Peaceman,  and  A.L.  Pozzi:  Numerical  calculations  of  multidimen¬ 
sional  miscible  displacement  by  the  method  of  characteristics,  Soc.  Pet.  Enq.  J.,  4, 
(1964),  26-36. 

[37]  P.  Hansbo  and  A.  Szepessy:  A  velocity-pressure  streamline  diffusion  finite  element 
method  for  the  incompressible  Navier-Stokes  equations,  Comp.  Meth.  Appl.  Mech.  En- 
grg.,  84,  (1990),  107-129. 

[38]  P.  Hansbo:  The  characteristic  streamline  diffusion  method  for  the  time-independent 
incompressible  Navier-Stokes  equations,  Comp.  Meth.  Appl.  Mech.  Engrg.,  99,  (1992), 
171-186. 

[39]  A.  Harten,  B.Engquist,  S.  Oslier,  and  S.  Chakravarthy:  Uniformly  high  order  accurate 
essentially  non-oscillatory  schemes,  III,  J.  Comput.  Phys.  71,  (1987),  231-241. 

[40]  R.W.  Healy  and  T.F.  Russell:  A  finite-volume  Eulerian-Lagrangian  localized  adjoint 
method  for  solution  of  the  advection-dispersion  equation,  Water  Resource  Research, 
29,  (1993),  2399-2413. 

[41]  R.W.  Healy  and  T.F.  Russell:  Solution  of  the  advection-dispersion  equation  in  two 
dimensions  by  a  finite-volume  Eulerian-Lagrangian  localized  adjoint  method,  Advances 
in  Water  Resources,  21,  (1998),  11-26. 

[42]  I.  Herrera,  R.E.  Ewing,  M.A.  Celia,  and  T.F.  Russell:  Eulerian-Lagrangian  localized 
adjoint  methods:  the  theoretical  framework,  Num.  Meth.  PDE’s,  9,  (1993),  431-458. 

[43]  J.M.  Hervouet:  Applications  of  the  method  of  characteristics  in  their  weak  formulation 
to  solving  two-dimensional  advection-equations  on  mesh  grids,  Computational  Tech¬ 
niques  for  Fluid  Flow ,  Recent  Advances  in  Numerical  Methods  in  Fluids ,  5,  Taylor  et 
al.  (eds.),  Pineidge  Press,  (1986),  149-185. 

[44]  T.J.R.  Hughes  and  A.N.  Brooks:  A  multidimensional  upwinding  scheme  with  no  cross- 
wind  diffusion,  Finite  Element  Methods  for  Convection  Dominated  Flows,  34,  Hughes 
(ed.),  ASME,  New  York,  (1979). 

[45]  T.J.R.  Hughes  and  M.  Mallet:  A  new  finite  element  formulation  for  computational  fluid 
dynamics:  III.  The  general  streamline  operator  for  multidimensional  advective-diffusive 
systems,  Comp.  Meth.  Appl.  Mech.  Engrg.,  58,  (1986),  305-328. 

[46]  T.J.R.  Hughes:  Multiscale  phenomena:  Green  functions,  the  Dirichlet-to-Neumann  for¬ 
mulation,  subgrid  scale  models,  bubbles,  and  the  origins  of  stabilized  methods.  Comp. 
Meth.  Appl.  Mech.  Engrg.,  127,  (1995),  387-401. 

[47]  C.  Johnson  and  U.  Navert:  An  analysis  of  some  finite  element  methods  for  advection- 
diffusion  equations,  Alexsson  et  al.  (eds.),  Analytical  and  Numerical  Approaches  to 
Asymptotic  Problems  in  Analysis ,  North-Holland,  (1981). 

[48]  C.  Johnson  and  J.  Pitkaranta:  An  analysis  of  discontinuous  Galerkin  methods  for  a 
scalar  hyperbolic  equation,  Math.  Comp.,  46,  (1986),  1-26. 

[49]  C.  Johnson:  Numerical  solutions  of  partial  differential  equations  by  the  finite  element 
method,  Cambridge  University  Press,  Cambridge,  1987. 

[50]  C.  Johnson,  A.  Szepessy,  and  P.  Hansbo:  On  the  convergence  of  shock-capturing  stream¬ 
line  diffusion  finite  element  methods  for  hyperbolic  conservation  laws,  Math.  Comp.,  54, 
(1990),  107-129. 


28 


[51]  C.  Johnson  and  A.  Szepessy:  Adaptive  finite  element  methods  for  conservation  laws 
based  on  a  posteriori  error  estimates,  Comm.  Pure  Appl.  Math.,  98,  (1995),  199-234. 

[52]  G.  Lube  and  D.  Weiss:  Stabilized  finite  element  methods  for  singularly  perturbed 
parabolic  problems,  Appl.  Num.  Math.,  17,  (1995),  431-459. 

[53]  N.  Lu:  A  semianalytical  method  of  path  line  computation  for  transient  finite-difference 
groundwater  flow  models,  Water  Resource  Research,  30,  (1994),  2449-2459. 

[54]  K.W.  Morton,  A.  Priestley,  and  E.  Stili:  Stability  of  the  Lagrangian-Galerkin  method 
with  nonexact  integration,  RAIRO  Ad2 AN,  22,  (1988),  123-151. 

[55]  S.P.  Neuman:  An  Eulerian-Lagrangian  numerical  scheme  for  the  dispersion-convection 
equation  using  conjugate  space-time  grids,  J.  Comp.  Phys.,  41,  (1981),  270-294. 

[56]  G.F.  Pinder  and  H.H.  Cooper:  A  numerical  technique  for  calculating  the  transient 
position  of  the  saltwater  front,  Water  Resource  Research ,  (1970),  875-882. 

[57]  O.  Pironneau:  On  the  transport-diffusion  algorithm  and  its  application  to  the  Navier- 
Stokes  equations,  Num.  Math.,  38,  (1982),  309-332. 

[58]  D.W.  Pollock:  Semianalytical  computation  of  path  lines  for  finite-difference  models, 
Ground  Water,  26,  (1988),  743-750. 

[59]  G.R.  Richter:  An  optimal-order  error  estimate  for  the  discontinuous  Galerkin  method, 
Math.  Comp.  50,  (1988),  75-88. 

[60]  T.F.  Russell:  Eulerian-Lagrangian  localized  adjoint  methods  for  advection-dominated 
problems,  Proceedings  of  the  13th  Dundee  Conference  on  Numerical  Analysis ,  Griffiths 
and  Watson  (eds.),  Pitmann  Research  Notes  in  Mathematics  Series,  228,  Longman 
Scientific  &  Technical,  Harlow,  U.  K.,  (1990),  206-228. 

[61]  T.F.  Russell,  and  R.V.  Trujillo:  Eulerian-Lagrangian  localized  adjoint  methods  with 
variable  coefficients  in  multiple  dimensions,  Gambolati,  et  al.  (ed.),  Computational 
Methods  in  Surface  Hydrology,  Springer- Verlag,  Berlin,  (1990),  357-363. 

[62]  A.L.  Scliafer-Perini  and  J.L.  Wilson:  Efficient  and  accurate  front  tracking  for  two- 
dimensional  groundwater  flow  models,  Water  Resource  Research,  27,  (1991),  1471-1485. 

[63]  C.  Shu  and  S.  Oslier:  Efficient  implementation  of  essentially  non-oscillatory  shock  cap¬ 
turing  schemes,  J.  Comp.  Phys.  77,  (1988),  439-471. 

[64]  J.E.  Vag,  H.  Wang,  and  H.K.  Dahle:  Eulerian-Lagrangian  localized  adjoint  methods 
for  systems  of  nonlinear  advective-diffusive-reactive  equations,  Advances  in  Water  Re¬ 
sources,  19,  (1996),  297-315. 

[65]  B.  van  Leer:  On  the  relation  between  the  upwind-differencing  schemes  of  Godunov, 
Engquist-Osher,  and  Roe,  SIAM  J.  Sci.  Stat.  Comp.  5,  (1984),  1-20. 

[66]  E.  Varoglu  and  W.D.L.  Finn:  Finite  elements  incorporating  characteristics  for  one¬ 
dimensional  diffusion-convection  equation,  J.  Comp.  Phys.,  34,  (1980),  371-389. 

[67]  H.  Wang:  Eulerian-Lagrangian  localized  adjoint  methods:  analyses,  numerical  imple¬ 
mentations  and  their  applications,  Ph.D.  Thesis,  Department  of  Mathematics,  Univer¬ 
sity  of  Wyoming,  1992. 

[68]  H.  Wang,  R.E.  Ewing,  and  T.F.  Russell:  ELLAM  for  variable-coefficient  convection- 
diffusion  problems  arising  in  groundwater  applications,  Computational  Methods  in  Wa¬ 
ter  Resources  IX,  1,  Computational  Mechanics  Publications  and  Elsevier  Applied  Sci¬ 
ence,  London  and  New  York,  (1992),  25-31. 


29 


[69]  H.  Wang  and  R.E.  Ewing:  Optimal-order  convergence  rates  for  ELLAM  for  reactive 
transport  and  contamination  in  groundwater,  Num.  Meth.  PDEs ,  11,  (1995),  1-31. 

[70]  H.  Wang,  R.E.  Ewing,  and  T.F.  Russell:  Eulerian-Lagrangian  localized  methods  for 
convection-diffusion  equal  ions  and  their  convergence  analysis,  IMA  J.  Num.  Anal. , 
15(3),  (1995a). 

[71]  H.  Wang,  R.E.  Ewing,  and  M.A.  Celia:  Eulerian-Lagrangian  localized  adjoint  method 
for  reactive  transport  with  biodegradation,  Num.  Meth.  PDEs,  11,  (1995b),  229-254. 

[72]  H.  Wang,  R.C.  Sharpley,  and  S.  Man:  An  ELLAM  scheme  for  advection-diffusion 
equations  in  multi-dimensions,  Aldama  et  al.  (eds.)  Computational  Methods  in  Wa¬ 
ter  Resources  XI,  2,  Computational  Mechanics  Publications,  Southampton  and  Boston, 
(1996),  99  -  106. 

[73]  H.  Wang,  M.A.  Al-Lawatia,  and  A.S.  Telyakovskiy:  Runge-Kutta  characteristic  methods 
for  first-order  linear  hyperbolic  equations,  to  appear  in  Num.  Meth.  PDEs. 

[74]  J.J.  Westerink  and  D.  Shea:  Consider  higher  degree  Petrov-Galerkin  methods  for  the 
solution  of  the  transient  convection-diffusion  equation,  Int.  J.  Num.  Meth.  Engrg.,  28, 
(1989),  1077-1101. 

[75]  M.F.  Wheeler  and  C.N.  Dawson:  An  operator-splitting  method  for  advection-diffusion- 
reaction  problems,  MAFELAP  Proceedings ,  6,  (Whiteman  ed.),  Academic  Press,  (1988), 
463-482. 

[76]  D.  Yang:  A  characteristic-mixed  method  with  dynamic  finite  element  space  for 
convection-dominated  diffusion  problems,  J.  Comp.  Appl.  Math.,  43,  (1992),  343-353. 

[77]  G.  Zhou:  An  adaptive  streamline  diffusion  finite  element  method  for  hyperbolic  systems 
in  gas  dynamics.  Ph.D.  Thesis,  Department  of  Mathematics,  University  of  Heidelberg, 
Heidelberg,  Germany,  1992. 

[78]  G.  Zhou:  A  local  L2-error  analysis  of  the  streamline  diffusion  method  for  nonstationary 
convection-diffusion  systems,  Mathematical  Modeling  and  Numerical  Analysis,  RAIRO 
25,  (1995),  577-603. 


30 


Analytical 


ELLAM-Euler 


ELLAM-RK 


BE-GAL 


SDM,  K=0.5 


K=0.01 


K=0.001 


Max. 

Min. 

CPU  (sec.) 
/time  step 

Overall  CPU 

0.8642 

0 

N/A 

N/A 

0.8302 

0 

17 

5  min.  39  sec. 

0.8630 

0 

10.2 

2  min.  33  sec. 

0.2546 

0 

15.4 

25  min.  40  sec. 

0.4517 

0 

13.85 

1  hr.  32 i  min. 

0.4523 

0 

31.75 

3  hr.  31 1  min. 

0.6006 

0 

13.65 

3  hr.  47 i  min. 

0.6015 

0 

30.85 

8  hr.  34^  min. 

6 

0.6964 

0 

13.6 

7  hr.  33-j  min. 

0.7401 

0 

13.55 

11  hr.  17^  min. 

0.2003 

-0.0002 

17.75 

29  min.  35  sec. 

0.3486 

-0.0003 

17.9 

1  hr.  59-j  min. 

0.3842 

0 

39.5 

4  hr.  23-j  min. 

0.4563 

-0.0008 

18.1 

5  hr.  l|  min. 

0.5035 

-0.0003 

40.1 

11  hr.  8-j-  min. 

0.5250 

-0.0015 

18.25 

10  hr.  8-|-  min. 

0.5554 

-0.0017 

18.45 

15  hr.  22  \  min. 

0.7861 

-0.1564 

15.05 

25  min.  5  sec. 

0.8438 

-0.0159 

13.95 

46  min. 

0.6197 

-0.0978 

17.95 

29  min.  55  sec. 

0.6412 

-0.0081 

18 

1  hr. 

unbounded 

unbounded 

N/A 

N/A 

0.8555 

-0.0002 

20 

1  hr.  6|  min. 

0.7089 

-0.0147 

96.9 

2  hr.  41.5  min. 

0.8250 

-0.0021 

96.9 

2  hr.  41.5  min. 

0.8281 

-0.0019 

96.9 

2  hr.  41.5  min. 

Figure  # 


5  (c)-(d) 


6  (a)-(b) 


6  (c)-(d) 


7  (a)-(b) 


8  (a) -(b) 


7  (c)-(d) 


8  (c)-(d) 


9  (a)-(b) 


10  (a)-(b) 


9  (c)-(d) 


10  (c)-(d) 


11  (a)-(b) 


11  (c)-(d) 


12  (a)-(b) 


12  (c)-(d) 


31 


(b)  Tost  functions:  for  GAL,  and  "  for  QPG  with  V  =  1,D  =  0.3  and  V  = 


(c)  Tost  functions:  for  GAL,  and  “ — “  for  CPG  with  Cu  =  0.5  and  Cu 

Fig  1.  Illustration  of  characteristic  tracking  and  test  functions 


